heatEquation_example.cpp

Annotated source file

Here is the full file examples/heatEquation_example.cpp. Clicking on a function or class name will lead you to its reference documentation.

#include <gismo.h>


using namespace gismo;


int main(int argc, char *argv[])
{
    bool plot = false;
    gsCmdLine cmd("Testing the heat equation.");
    cmd.addSwitch("plot", "Plot the result in ParaView.", plot);
    try { cmd.getValues(argc,argv); } catch (int rv) { return rv; }

    // Source function
    gsConstantFunction<> f(0,2);
    gsInfo<<"Source function is: "<< f << "\n";

    // Define Geometry, must be a gsMultiPatch object
    gsMultiPatch<> patches(*gsNurbsCreator<>::BSplineSquareDeg(2));
    patches.computeTopology();

    // Boundary conditions
    gsBoundaryConditions<> bcInfo;
    gsConstantFunction<> g_N(1,2); // Neumann
    gsConstantFunction<> g_D(0,2); // (Dirichlet
    bcInfo.addCondition(0, boundary::west,  condition_type::neumann  , &g_N);
    bcInfo.addCondition(0, boundary::east,  condition_type::dirichlet, &g_D);
    bcInfo.addCondition(0, boundary::north, condition_type::dirichlet, &g_D);
    bcInfo.addCondition(0, boundary::south, condition_type::dirichlet, &g_D);

    gsMultiBasis<> refine_bases( patches );
    // Number for h-refinement of the computational (trial/test) basis.
    int numRefine  = 2;

    // Number for p-refinement of the computational (trial/test) basis.
    int numElevate = 0;

    // Elevate and p-refine the basis to order k + numElevate
    // where k is the highest degree in the bases
    if ( numElevate > -1 )
    {
        // Find maximum degree with respect to all the variables
        int tmp = refine_bases.maxDegree(0);
        for (short_t j = 1; j < patches.parDim(); ++j )
            if ( tmp < refine_bases.maxDegree(j) )
                tmp = refine_bases.maxDegree(j);

        // Elevate all degrees uniformly
        tmp += numElevate;
        refine_bases.setDegree(tmp);
    }

    // h-refine the basis
    for (int i = 0; i < numRefine; ++i)
        refine_bases.uniformRefine();

    // Determines the theta-scheme used for time integration
    // (eg. Forward/backward Euler or Crank Nicolson(theta=0.5)
    real_t theta = 0.5;

    gsPoissonPde<> pde(patches, bcInfo, f);
    // Assembler (constructs stationary matrix and right-hand side vector)
    gsPoissonAssembler<> stationary(pde, refine_bases);
    stationary.options().setInt("DirichletStrategy", dirichlet::elimination);
    stationary.options().setInt("InterfaceStrategy", iFace::glue);
    gsHeatEquation<real_t> assembler(stationary);
    assembler.setTheta(theta);
    gsInfo<<assembler.options()<<"\n";

    // A Conjugate Gradient linear solver with a diagonal (Jacobi) preconditionner
    gsSparseSolver<>::CGDiagonal solver;

    // Generate system matrix and load vector
    gsInfo<<"Assembling mass and stiffness...\n";
    assembler.assemble();

    gsMatrix<> Sol, Rhs;
    int ndof = assembler.numDofs();
    real_t endTime = 0.1;
    int numSteps = 40;
    Sol.setZero(ndof, 1); // Initial solution

    real_t Dt = endTime / numSteps ;

    const std::string baseName("heat_eq_solution");
    gsParaviewCollection collection(baseName);
    collection.options().setInt("numPoints", 1000);
    collection.options().setInt("precision", 5);

    if ( plot )
    {
        //sol = assembler.constructSolution(Sol); // same as next line
        gsField<> sol = stationary.constructSolution(Sol);

        collection.newTimeStep(&patches);
        collection.addField(sol, "Temperature");
        collection.saveTimeStep();
    }

    for ( int i = 1; i<=numSteps; ++i) // for all timesteps
    {
        // Compute the system for the timestep i (rhs is assumed constant wrt time)
        assembler.nextTimeStep(Sol, Dt);
        gsInfo<<"Solving timestep "<< i*Dt<<".\n";

        // Solve for current timestep, overwrite previous solution
        Sol = solver.compute( assembler.matrix() ).solve( assembler.rhs() );

        // Obtain current solution as an isogeometric field
        //sol = assembler.constructSolution(Sol); // same as next line
        gsField<> sol = stationary.constructSolution(Sol);

        if ( plot )
        {
            // Plot the snapshot to paraview
            collection.newTimeStep(&patches);
            collection.addField(sol, "Temperature");
            collection.saveTimeStep();
        }
    }

    //gsInfo<< " time = "<<endTime<<"\n";

    if ( plot )
    {
        collection.save();
        gsFileManager::open("heat_eq_solution.pvd");
    }
    else
        gsInfo << "Done. No output created, re-run with --plot to get a ParaView "
                  "file containing the solution.\n";

    return  EXIT_SUCCESS;
}