ex2.cpp

Go to the documentation of this file.

#include <iostream>

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

#include <adaptiveTriangMesh.h>
#include <configuratorsShellFE.h>
#include <EigenDataContainer.h>
#include <quadratureShellFE.h>
#include <legacyVtkWriter.h>
#include <shellHandler.h>
#include <triangleShellFE.h>
#include "errorEstimator.h"

// typedefs for better code readability
typedef double RealType;
typedef shellFE::ShellElementWithTangentSpaceAtVertex<DataTypeContainerShellFE> TriangleType;
typedef shellFE::AdaptiveTriangMesh<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: ex2 meshFile resultFile [numRefinements]\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";
    cout << "\tnumRefinements\tNumber of times the mesh is refined"
            << " (default: 10)\n\n";
}

int main(int argc, char **argv)
{
    unsigned int numRefinements = 10;
    switch (argc) {
    case 3: {
        // First two arguments are handled below
        break;
    }
    case 4: {
        numRefinements = atoi(argv[3]);
        break;
    }
    default: {
        printHelp();
        return 1;
    }
    }

    // Load mesh
    cout << "* Loading grid" << endl;
    MeshType mesh(argv[1]);

    cout << "* Refining " << numRefinements << " times ";
    cout << "  [";
    cout.flush();
    for (unsigned int r = 0; r < numRefinements; ++r) {
        // Estimate error and refine mesh
        ErrorEstimator<ConfiguratorType> errorEstimator(mesh);
        errorEstimator.markAndRefineTriangles();

        // Save the mesh and solution as legacy VTK
        shellFE::LegacyVtkWriter<MeshType> vtkWriter(mesh);
        std::stringstream outputFileName;
        outputFileName << argv[2] << r << ".vtk";
        vtkWriter.save(outputFileName.str().c_str());

        cout << '|';
        cout.flush();
    }
    cout << "]\n\n";

    return 0;
}