biharmonic2_example.cpp

The boundary value problem for the first biharmonic equation is defined as follows:

\begin{eqnarray*} \Delta^2 u &=& f \quad \text{in} \quad \Omega, \\ u &=& g_1 \quad \text{on} \quad \partial \Omega, \\ \mathbf{n} \cdot \nabla u &=& g_2 \quad \text{on} \quad \partial \Omega, \end{eqnarray*}

where \(\mathbf{n}\) is the outward pointing unit normal vector. The boundary value problem for the second biharmonic equation is defined as follows:

\begin{eqnarray*} \Delta^2 u &=& f \quad \text{in} \quad \Omega, \\ u &=& g_1 \quad \text{on} \quad \partial \Omega, \\ \Delta u &=& g_2 \quad \text{on} \quad \partial \Omega. \end{eqnarray*}

The system to assemble a system of equations based on this weak formulation is constructed in biharmonic2_example.cpp using the gsExprAssembler. The biharmonic problem is only solvable on domains with a single patch.

Then we define our command line options. For example, we use the option -f or –file to set the path to the file that contains our geometry. If we wish to use the second biharmonic problem, then we use the option –second to run the biharmonic problem with \(u=g_1\) and \(\Delta u=g_2\)

    bool plot = false;

    index_t numRefine  = 3;
    index_t degree = 3;
    index_t smoothness = 2;
    bool last = false;
    bool second = false;

    std::string fn;

    gsCmdLine cmd("Example for solving the biharmonic problem (single patch only).");
    cmd.addInt("p", "degree","Set discrete polynomial degree", degree);
    cmd.addInt("s", "smoothness", "Set discrete regularity",  smoothness);
    cmd.addInt("l", "refinementLoop", "Number of refinement steps",  numRefine);
    cmd.addString("f", "file", "Input geometry file (with .xml)", fn);

    cmd.addSwitch("last", "Solve problem only on the last level of h-refinement", last);
    cmd.addSwitch("plot", "Create a ParaView visualization file with the solution", plot);

    cmd.addSwitch("second", "Compute second biharmonic problem with u = g1 and Delta u = g2 "
                            "(default first biharmonic problem: u = g1 and partial_n u = g2)", second);
    try { cmd.getValues(argc,argv); } catch (int rv) { return rv; }

Then we read the geometry: Either it can be from a file that contains the geometry or we use the default one.

    gsMultiPatch<> mp;
    if (fn.empty())
        mp = gsMultiPatch<>( *gsNurbsCreator<>::BSplineFatQuarterAnnulus() );
    else
    {
        gsInfo << "Filedata: " << fn << "\n";
        gsReadFile<>(fn, mp);
    }
    mp.clearTopology();
    mp.computeTopology();

    if (mp.nPatches() != 1)
    {
        gsInfo << "The geometry has more than one patch. Run the code with a single patch!\n";
        return EXIT_FAILURE;
    }

Before we run the problem, we set the degree and/or refine the basis. This is done in the following snippet:

    gsMultiBasis<> basis(mp, false);//true: poly-splines (not NURBS)

    // Elevate and p-refine the basis to order p + numElevate
    // where p is the highest degree in the bases
    basis.setDegree(degree); // preserve smoothness

    // h-refine each basis
    if (last)
    {
        for (index_t r =0; r < numRefine; ++r)
            basis.uniformRefine(1, degree-smoothness);
        numRefine = 0;
    }

Then we set the boundary conditions. Depending on which biharmonic problem we choose, each boundary will be is set to the corresponding functions.

    gsBoundaryConditions<> bc;
    for (gsMultiPatch<>::const_biterator bit = mp.bBegin(); bit != mp.bEnd(); ++bit)
    {
        // Laplace
        gsFunctionExpr<> laplace ("-16*pi*pi*(2*cos(4*pi*x)*cos(4*pi*y) - cos(4*pi*x) - cos(4*pi*y))",2);

        // Neumann
        gsFunctionExpr<> sol1der("-4*pi*(cos(4*pi*y) - 1)*sin(4*pi*x)",
                                 "-4*pi*(cos(4*pi*x) - 1)*sin(4*pi*y)", 2);

        bc.addCondition(*bit, condition_type::dirichlet, ms);
        if (second)
            bc.addCondition(*bit, condition_type::laplace, laplace);
        else
            bc.addCondition(*bit, condition_type::neumann, sol1der);
    }
    bc.setGeoMap(mp);
    gsInfo << "Boundary conditions:\n" << bc << "\n";

Now it is time to look at the setup of our boundary value problem. In the code snippet below, we first define our gsExprAssembler. The inputs (1,1) here mean that we have 1 test function space and 1 trial function space. Furthermore, few typedef s can be seen. These are just there to prevent lengthy types. The assembler is fed with the basis using the gsExprAssembler::setIntegrationElements() function. When the basis is defined, a gsExprEvaluator is constructed for later use. The geometry is connected to the gsExprAssembler via the gsExprAssembler::getMap() function and we create an expression for the space via gsExprAssembler::getSpace(). Functions for the right-hand side and the manufactured solution are constructed via gsExprAssembler::getCoef() and gsExprEvaluator::getVariable() functions, respectively. Note that the manufactured solution is defined in the gsExprEvaluator because it will not be used in assembly of the system. The problem setup finishes by defining the solution object, which is constructed via gsExprAssembler::getSolution() which requires a space and a reference to a vector that will contain the solution coefficients later.

    gsExprAssembler<real_t> A(1,1);
    //gsInfo<<"Active options:\n"<< A.options() <<"\n";

    // Elements used for numerical integration
    A.setIntegrationElements(basis);
    gsExprEvaluator<real_t> ev(A);

    // Set the geometry map
    auto G = A.getMap(mp);

    // Set the source term
    auto ff = A.getCoeff(f, G); // Laplace example

    // Set the discretization space
    auto u = A.getSpace(basis);

    // Solution vector and solution variable
    gsMatrix<real_t> solVector;
    auto u_sol = A.getSolution(u, solVector);

    // Recover manufactured solution
    auto u_ex = ev.getVariable(ms, G);

The next step is to actually solve the system at hand. In this particular example, this is done within a loop of uniform refinement (unless the last option is selected) in order to study convergence behaviour of the method.

As seen in the snippet below, we first define the solver that we will use (based on the Eigen-library) and we define the vectors to store the errors.

Then, within the loop, we first refine the basis uniformly.

For the biharmonic problem, we fix the dofs in the function setMapperForBiharmonic(), which is defined in the beginning of the file (see the full file at the end).

Then the gsExprAssembler<>::space::setupMapper() function is used. This function initialize the matrix and the rhs depending on the mapper.

Then we compute the eliminated boundary values with the L2-projection where we call setDirichletNeumannValuesL2Projection(). Again this function is defined at the beginning of the file.

Then the gsExprAssembler is initialized using gsExprAssembler::initSystem() and using gsExprAssembler::assemble() the matrix expression (first argument) and the rhs expression (second argument) are assembled to a system.

When the System is assembled, the system is factorized using gsSparseSolver::compute() and solved via gsSparseSolver::solve(). The errors are then computed via the gsExprEvaluator, where it should be noted that the solution object is appearing, which uses the solution vector and the space as defined earlier.

    gsSparseSolver<real_t>::SimplicialLDLT solver;

    gsVector<real_t> l2err(numRefine+1), h1err(numRefine+1), h2err(numRefine+1),
            dofs(numRefine+1), meshsize(numRefine+1);
    gsInfo<< "(dot1=assembled, dot2=solved, dot3=got_error)\n"
             "\nDoFs: ";
    double setup_time(0), ma_time(0), slv_time(0), err_time(0);
    gsStopwatch timer;
    for (index_t r=0; r<=numRefine; ++r)
    {
        // Refine uniform once
        basis.uniformRefine(1,degree -smoothness);
        meshsize[r] = basis.basis(0).getMaxCellLength();

        // Setup the Mapper
        gsDofMapper map;
        setMapperForBiharmonic(bc, basis,map);

        // Setup the system
        u.setupMapper(map);
        setDirichletNeumannValuesL2Projection(mp, basis, bc, u);

        // Initialize the system
        A.initSystem();
        setup_time += timer.stop();

        gsInfo<< A.numDofs() <<std::flush;

        timer.restart();
        // Compute the system matrix and right-hand side
        A.assemble(ilapl(u, G) * ilapl(u, G).tr() * meas(G),u * ff * meas(G));

        // Enforce Laplace conditions to right-hand side
        auto g_L = A.getBdrFunction(G); // Set the laplace bdy value
        //auto g_L = A.getCoeff(laplace, G);
        A.assembleBdr(bc.get("Laplace"), (igrad(u, G) * nv(G)) * g_L.tr() );

        dofs[r] = A.numDofs();
        ma_time += timer.stop();
        gsInfo << "." << std::flush;// Assemblying done

        timer.restart();
        solver.compute( A.matrix() );
        solVector = solver.solve(A.rhs());

        slv_time += timer.stop();
        gsInfo << "." << std::flush; // Linear solving done

        timer.restart();
        //linferr[r] = ev.max( f-s ) / ev.max(f);

        l2err[r]= math::sqrt( ev.integral( (u_ex - u_sol).sqNorm() * meas(G) ) ); // / ev.integral(ff.sqNorm()*meas(G)) );
        h1err[r]= l2err[r] +
                  math::sqrt(ev.integral( ( igrad(u_ex) - igrad(u_sol,G) ).sqNorm() * meas(G) )); // /ev.integral( igrad(ff).sqNorm()*meas(G) ) );

        h2err[r]= h1err[r] +
                  math::sqrt(ev.integral( ( ihess(u_ex) - ihess(u_sol,G) ).sqNorm() * meas(G) )); // /ev.integral( ihess(ff).sqNorm()*meas(G) )

        err_time += timer.stop();
        gsInfo << ". " << std::flush; // Error computations done
    } //for loop

When the loop is runnig several times, we can compute the L2-error, H1-error and H2-error. The expected rates should be p+1 for the L2-error, p for the H1-error and p-1 for the H2-error.

    gsInfo << "\nL2 error: "<<std::scientific<<std::setprecision(3)<<l2err.transpose()<<"\n";
    gsInfo << "H1 error: "<<std::scientific<<h1err.transpose()<<"\n";
    gsInfo << "H2 error: "<<std::scientific<<h2err.transpose()<<"\n";

    if (numRefine>0)
    {
        gsInfo<< "\nEoC (L2): " << std::fixed<<std::setprecision(2)
              <<  ( l2err.head(numRefine).array()  /
                    l2err.tail(numRefine).array() ).log().transpose() / std::log(2.0)
              <<"\n";

        gsInfo<<   "EoC (H1): "<< std::fixed<<std::setprecision(2)
              <<( h1err.head(numRefine).array() /
                  h1err.tail(numRefine).array() ).log().transpose() / std::log(2.0) <<"\n";

        gsInfo<<   "EoC (H2): "<< std::fixed<<std::setprecision(2)
              <<( h2err.head(numRefine).array() /
                  h2err.tail(numRefine).array() ).log().transpose() / std::log(2.0) <<"\n";
    }

If we set the flag –plot then the export to ParaView is done in the snipped below. The output file is solution.pvd.

    if (plot)
    {
        // Write approximate and exact solution to paraview files
        gsInfo << "Plotting in Paraview...\n";
        ev.options().setSwitch("plot.elements", true);
        ev.writeParaview( u_sol   , G, "solution");
        gsInfo << "Saved with solution.pvd \n";
    }
    else
        gsInfo << "Done. No output created, re-run with --plot to get a ParaView "
                  "file containing the solution.\n";

Annotated source file

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

# include <gismo.h>

using namespace gismo;


void setMapperForBiharmonic(gsBoundaryConditions<> & bc, gsMultiBasis<> & basis, gsDofMapper & mapper)
{
    mapper.init(basis);

    for (gsBoxTopology::const_iiterator it = basis.topology().iBegin();
         it != basis.topology().iEnd(); ++it) // C^0 at the interface
    {
        basis.matchInterface(*it, mapper);
    }

    gsMatrix<index_t> bnd;
    for (typename gsBoundaryConditions<real_t>::const_iterator
                 it = bc.begin("Dirichlet"); it != bc.end("Dirichlet"); ++it)
    {
        bnd = basis.basis(it->ps.patch).boundary(it->ps.side());
        mapper.markBoundary(it->ps.patch, bnd, 0);
    }

    for (typename gsBoundaryConditions<real_t>::const_iterator
                 it = bc.begin("Neumann"); it != bc.end("Neumann"); ++it)
    {
        bnd = basis.basis(it->ps.patch).boundaryOffset(it->ps.side(),1);
        mapper.markBoundary(it->ps.patch, bnd, 0);
    }
    mapper.finalize();
}

void setDirichletNeumannValuesL2Projection(gsMultiPatch<> & mp, gsMultiBasis<> & basis, gsBoundaryConditions<> & bc, const expr::gsFeSpace<real_t> & u)
{
    const gsDofMapper & mapper = u.mapper();
    gsDofMapper mapperBdy(basis, u.dim());
    for (gsBoxTopology::const_iiterator it = basis.topology().iBegin();
         it != basis.topology().iEnd(); ++it) // C^0 at the interface
    {
        basis.matchInterface(*it, mapperBdy);
    }
    for (size_t np = 0; np < mp.nPatches(); np++)
    {
        gsMatrix<index_t> bnd = mapper.findFree(np);
        mapperBdy.markBoundary(np, bnd, 0);
    }
    mapperBdy.finalize();

    gsExprAssembler<real_t> A(1,1);
    A.setIntegrationElements(basis);

    auto G = A.getMap(mp);
    auto uu = A.getSpace(basis);
    auto g_bdy = A.getBdrFunction(G);

    uu.setupMapper(mapperBdy);
    gsMatrix<real_t> & fixedDofs_A = const_cast<expr::gsFeSpace<real_t>&>(uu).fixedPart();
    fixedDofs_A.setZero( uu.mapper().boundarySize(), 1 );

    real_t lambda = 1e-5;

    A.initSystem();
    A.assembleBdr(bc.get("Dirichlet"), uu * uu.tr() * meas(G));
    A.assembleBdr(bc.get("Dirichlet"), uu * g_bdy * meas(G));
    A.assembleBdr(bc.get("Neumann"),
                  lambda * (igrad(uu, G) * nv(G).normalized()) * (igrad(uu, G) * nv(G).normalized()).tr() * meas(G));
    A.assembleBdr(bc.get("Neumann"),
                  lambda *  (igrad(uu, G) * nv(G).normalized()) * (g_bdy.tr() * nv(G).normalized()) * meas(G));

    gsSparseSolver<real_t>::SimplicialLDLT solver;
    solver.compute( A.matrix() );
    gsMatrix<real_t> & fixedDofs = const_cast<expr::gsFeSpace<real_t>& >(u).fixedPart();
    gsMatrix<real_t> fixedDofs_temp = solver.solve(A.rhs());

    // Reordering the dofs of the boundary
    fixedDofs.setZero(mapper.boundarySize(),1);
    index_t sz = 0;
    for (size_t np = 0; np < mp.nPatches(); np++)
    {
        gsMatrix<index_t> bnd = mapperBdy.findFree(np);
        bnd.array() += sz;
        for (index_t i = 0; i < bnd.rows(); i++)
        {
            index_t ii = mapperBdy.asVector()(bnd(i,0));
            fixedDofs(mapper.global_to_bindex(mapper.asVector()(bnd(i,0))),0) = fixedDofs_temp(ii,0);
        }
        sz += mapperBdy.patchSize(np,0);
    }
}


int main(int argc, char *argv[])
{
    bool plot = false;

    index_t numRefine  = 3;
    index_t degree = 3;
    index_t smoothness = 2;
    bool last = false;
    bool second = false;

    std::string fn;

    gsCmdLine cmd("Example for solving the biharmonic problem (single patch only).");
    cmd.addInt("p", "degree","Set discrete polynomial degree", degree);
    cmd.addInt("s", "smoothness", "Set discrete regularity",  smoothness);
    cmd.addInt("l", "refinementLoop", "Number of refinement steps",  numRefine);
    cmd.addString("f", "file", "Input geometry file (with .xml)", fn);

    cmd.addSwitch("last", "Solve problem only on the last level of h-refinement", last);
    cmd.addSwitch("plot", "Create a ParaView visualization file with the solution", plot);

    cmd.addSwitch("second", "Compute second biharmonic problem with u = g1 and Delta u = g2 "
                            "(default first biharmonic problem: u = g1 and partial_n u = g2)", second);
    try { cmd.getValues(argc,argv); } catch (int rv) { return rv; }

    gsMultiPatch<> mp;
    if (fn.empty())
        mp = gsMultiPatch<>( *gsNurbsCreator<>::BSplineFatQuarterAnnulus() );
    else
    {
        gsInfo << "Filedata: " << fn << "\n";
        gsReadFile<>(fn, mp);
    }
    mp.clearTopology();
    mp.computeTopology();

    if (mp.nPatches() != 1)
    {
        gsInfo << "The geometry has more than one patch. Run the code with a single patch!\n";
        return EXIT_FAILURE;
    }

    gsFunctionExpr<>f("256*pi*pi*pi*pi*(4*cos(4*pi*x)*cos(4*pi*y) - cos(4*pi*x) - cos(4*pi*y))",2);
    gsInfo << "Source function: " << f << "\n";

    gsFunctionExpr<> ms("(cos(4*pi*x) - 1) * (cos(4*pi*y) - 1)",2);
    gsInfo << "Exact function: " << ms << "\n";

    gsMultiBasis<> basis(mp, false);//true: poly-splines (not NURBS)

    // Elevate and p-refine the basis to order p + numElevate
    // where p is the highest degree in the bases
    basis.setDegree(degree); // preserve smoothness

    // h-refine each basis
    if (last)
    {
        for (index_t r =0; r < numRefine; ++r)
            basis.uniformRefine(1, degree-smoothness);
        numRefine = 0;
    }

    gsBoundaryConditions<> bc;
    for (gsMultiPatch<>::const_biterator bit = mp.bBegin(); bit != mp.bEnd(); ++bit)
    {
        // Laplace
        gsFunctionExpr<> laplace ("-16*pi*pi*(2*cos(4*pi*x)*cos(4*pi*y) - cos(4*pi*x) - cos(4*pi*y))",2);

        // Neumann
        gsFunctionExpr<> sol1der("-4*pi*(cos(4*pi*y) - 1)*sin(4*pi*x)",
                                 "-4*pi*(cos(4*pi*x) - 1)*sin(4*pi*y)", 2);

        bc.addCondition(*bit, condition_type::dirichlet, ms);
        if (second)
            bc.addCondition(*bit, condition_type::laplace, laplace);
        else
            bc.addCondition(*bit, condition_type::neumann, sol1der);
    }
    bc.setGeoMap(mp);
    gsInfo << "Boundary conditions:\n" << bc << "\n";

    gsExprAssembler<real_t> A(1,1);
    //gsInfo<<"Active options:\n"<< A.options() <<"\n";

    // Elements used for numerical integration
    A.setIntegrationElements(basis);
    gsExprEvaluator<real_t> ev(A);

    // Set the geometry map
    auto G = A.getMap(mp);

    // Set the source term
    auto ff = A.getCoeff(f, G); // Laplace example

    // Set the discretization space
    auto u = A.getSpace(basis);

    // Solution vector and solution variable
    gsMatrix<real_t> solVector;
    auto u_sol = A.getSolution(u, solVector);

    // Recover manufactured solution
    auto u_ex = ev.getVariable(ms, G);

#ifdef _OPENMP
    gsInfo << "Available threads: "<< omp_get_max_threads() <<"\n";
#endif

    gsSparseSolver<real_t>::SimplicialLDLT solver;

    gsVector<real_t> l2err(numRefine+1), h1err(numRefine+1), h2err(numRefine+1),
            dofs(numRefine+1), meshsize(numRefine+1);
    gsInfo<< "(dot1=assembled, dot2=solved, dot3=got_error)\n"
             "\nDoFs: ";
    double setup_time(0), ma_time(0), slv_time(0), err_time(0);
    gsStopwatch timer;
    for (index_t r=0; r<=numRefine; ++r)
    {
        // Refine uniform once
        basis.uniformRefine(1,degree -smoothness);
        meshsize[r] = basis.basis(0).getMaxCellLength();

        // Setup the Mapper
        gsDofMapper map;
        setMapperForBiharmonic(bc, basis,map);

        // Setup the system
        u.setupMapper(map);
        setDirichletNeumannValuesL2Projection(mp, basis, bc, u);

        // Initialize the system
        A.initSystem();
        setup_time += timer.stop();

        gsInfo<< A.numDofs() <<std::flush;

        timer.restart();
        // Compute the system matrix and right-hand side
        A.assemble(ilapl(u, G) * ilapl(u, G).tr() * meas(G),u * ff * meas(G));

        // Enforce Laplace conditions to right-hand side
        auto g_L = A.getBdrFunction(G); // Set the laplace bdy value
        //auto g_L = A.getCoeff(laplace, G);
        A.assembleBdr(bc.get("Laplace"), (igrad(u, G) * nv(G)) * g_L.tr() );

        dofs[r] = A.numDofs();
        ma_time += timer.stop();
        gsInfo << "." << std::flush;// Assemblying done

        timer.restart();
        solver.compute( A.matrix() );
        solVector = solver.solve(A.rhs());

        slv_time += timer.stop();
        gsInfo << "." << std::flush; // Linear solving done

        timer.restart();
        //linferr[r] = ev.max( f-s ) / ev.max(f);

        l2err[r]= math::sqrt( ev.integral( (u_ex - u_sol).sqNorm() * meas(G) ) ); // / ev.integral(ff.sqNorm()*meas(G)) );
        h1err[r]= l2err[r] +
                  math::sqrt(ev.integral( ( igrad(u_ex) - igrad(u_sol,G) ).sqNorm() * meas(G) )); // /ev.integral( igrad(ff).sqNorm()*meas(G) ) );

        h2err[r]= h1err[r] +
                  math::sqrt(ev.integral( ( ihess(u_ex) - ihess(u_sol,G) ).sqNorm() * meas(G) )); // /ev.integral( ihess(ff).sqNorm()*meas(G) )

        err_time += timer.stop();
        gsInfo << ". " << std::flush; // Error computations done
    } //for loop

    timer.stop();
    gsInfo << "\n\nTotal time: "<< setup_time+ma_time+slv_time+err_time <<"\n";
    gsInfo << "     Setup: "<< setup_time <<"\n";
    gsInfo << "  Assembly: "<< ma_time    <<"\n";
    gsInfo << "   Solving: "<< slv_time   <<"\n";
    gsInfo << "     Norms: "<< err_time   <<"\n";
    gsInfo << " Mesh-size: "<< meshsize.transpose() << "\n";
    gsInfo << "      Dofs: "<<dofs.transpose() << "\n";

    gsInfo << "\nL2 error: "<<std::scientific<<std::setprecision(3)<<l2err.transpose()<<"\n";
    gsInfo << "H1 error: "<<std::scientific<<h1err.transpose()<<"\n";
    gsInfo << "H2 error: "<<std::scientific<<h2err.transpose()<<"\n";

    if (numRefine>0)
    {
        gsInfo<< "\nEoC (L2): " << std::fixed<<std::setprecision(2)
              <<  ( l2err.head(numRefine).array()  /
                    l2err.tail(numRefine).array() ).log().transpose() / std::log(2.0)
              <<"\n";

        gsInfo<<   "EoC (H1): "<< std::fixed<<std::setprecision(2)
              <<( h1err.head(numRefine).array() /
                  h1err.tail(numRefine).array() ).log().transpose() / std::log(2.0) <<"\n";

        gsInfo<<   "EoC (H2): "<< std::fixed<<std::setprecision(2)
              <<( h2err.head(numRefine).array() /
                  h2err.tail(numRefine).array() ).log().transpose() / std::log(2.0) <<"\n";
    }

    if (plot)
    {
        // Write approximate and exact solution to paraview files
        gsInfo << "Plotting in Paraview...\n";
        ev.options().setSwitch("plot.elements", true);
        ev.writeParaview( u_sol   , G, "solution");
        gsInfo << "Saved with solution.pvd \n";
    }
    else
        gsInfo << "Done. No output created, re-run with --plot to get a ParaView "
                  "file containing the solution.\n";

    return  EXIT_SUCCESS;
}