ex1.cpp

Go to the documentation of this file.

#include <iostream>

#include <Eigen/Dense>
#include <Eigen/Sparse>
#include <Eigen/SparseCholesky>

#include <configuratorsShellFE.h>
#include <EigenDataContainer.h>
#include <quadratureShellFE.h>
#include <legacyVtkWriter.h>
#include <shellHandler.h>
#include <triangMesh.h>
#include <triangleShellFE.h>
#include "rhs.h"
#include "stiffnessMatrixIntegrator.h"

// typedefs for better code readability
typedef double RealType;
typedef shellFE::ShellElementWithTangentSpaceAtVertex<DataTypeContainerShellFE> TriangleType;
typedef shellFE::TriangMesh<DataTypeContainerShellFE, TriangleType >            MeshType;
typedef shellFE::CenterOfMassQuadrature<RealType, typename DataTypeContainerShellFE::DomVecType> QuadType;
typedef shellFE::UnitTriangMeshConfiguratorP1<DataTypeContainerShellFE, MeshType, QuadType> ConfiguratorType;
typedef typename ConfiguratorType::VectorType VectorType;
typedef typename ConfiguratorType::SparseMatrixType SparseMatrixType;

void printHelp()
{
    cout << "Usage: ex1 meshFile resultFile\n";
    cout << "\tmeshFile\tLegacy VTK file defining the Finite Element mesh "
            << "to use\n";
    cout << "\tresultFile\tLegacy VTK file to write the result to\n\n";
}

int main(int argc, char **argv)
{
    if (argc != 3) {
        printHelp();
        return 1;
    }

    // Load a mesh
    MeshType mesh(argv[1]);
    // set neighbor relations
    mesh.makeNeighbour();

    // The configurator handles information about basis functions,
    // quadrature points etc.
    ConfiguratorType conf(mesh);

    // The shellHandler handles additional information about the mesh
    // including neighbor relations and the mapping from the unit triangle
    // to the elements of the mesh
    shellFE::ShellHandler<ConfiguratorType> shellHandler(conf);

    // Construct the right hand side
    VectorType rhs(conf.getNumGlobalDofs());
    createMaskedConstantRHS<ConfiguratorType>(conf, shellHandler, 1.0, rhs, true);

    // Construct the system matrix
    SparseMatrixType systemMat(conf.getNumGlobalDofs(), conf.getNumGlobalDofs());
    StiffnessMatrixIntegrator<ConfiguratorType> integrand(conf);
    // Assemble the system matrix
    integrand.assembleDirichlet<SparseMatrixType> (systemMat, shellHandler.getDirichletMask());

    // Construct a linear solver (Cholesky solver)
    Eigen::SimplicialLLT<SparseMatrixType> solver;
    // Decompose the matrix
    solver.compute(systemMat);
    // Solve the linear system
    VectorType sol = solver.solve(rhs);

    // Save the mesh and solution as legacy VTK
    LegacyVtkWriter<MeshType> vtkWriter ( mesh );
    vtkWriter.addScalarData(sol, "solution", VERTEX_DATA);
    vtkWriter.addScalarData(rhs, "rhs", VERTEX_DATA);
    vtkWriter.save(argv[2]);

    return 0;
}