ieti_example.cpp

We consider the Poisson boundary value problem (for simplicity here only with Dirichlet conditions):

\begin{eqnarray*} -\Delta u &=& f \quad \text{in} \quad \Omega, \\ u &=& g \quad \text{on} \quad \partial\Omega. \end{eqnarray*}

Its variational form is: Find \( u\in V_g \) such that

\[ \int_\Omega \nabla u \nabla v d x = \int_\Omega f v dx \]

for all \( v\in V_0 \), which corresponds to the linear system

\[ A \underline u = \underline f. \]

The computational domain is a non-overlapping multi-patch domain:

\[ \overline{ \Omega } = \bigcup_{k=1}^K \overline{ \Omega_k } \]

\[ \Omega_k \cap \Omega_l = \emptyset \quad \text{ for all } k\not=l. \]

On each patch, we have a isogeometric function space \( V_k \). We assume that the function spaces are fully matching. The global function space is

\[ V_h = \{ v\in V \;:\; v|_{\Omega_k} \in V_k \} \]

where \( V \) is the proper continuous space.

After reading the command line arguments, we read in the geometry file and modify the geometry as desired:

    gsMultiPatch<>::uPtr mpPtr = gsReadFile<>(geometry);

    gsMultiPatch<>& mp = *mpPtr;
 
    for (index_t i=0; i<splitPatches; ++i)
    {
        gsInfo << "split patches uniformly... " << std::flush;
        mp = mp.uniformSplit();
    }

Then, we set up the right-hand-side f and a gsBoundaryConditions object based on the command line arguments.

    // Right-hand-side
    gsFunctionExpr<> f( "2*sin(x)*cos(y)", mp.geoDim() );
 
    // Dirichlet function
    gsFunctionExpr<> gD( "sin(x)*cos(y)", mp.geoDim() );
 
    // Neumann
    gsConstantFunction<> gN( 1.0, mp.geoDim() );
 
    gsBoundaryConditions<> bc;

Then, we extract the basis from the multi patch object to obtain a multi basis.

    gsMultiBasis<> mb(mp);

We set the spline degree and the grid size as desired.

    for ( size_t i = 0; i < mb.nBases(); ++ i )
        mb[i].setDegreePreservingMultiplicity(degree);
 
    for ( index_t i = 0; i < refinements; ++i )
        mb.uniformRefine();

Now, we can discuss the IETI system.

We are considering patch-wise formulations of the variational formulation:

\[ \int_{\Omega_k} \nabla u \nabla v d x = \int_{\Omega_k} f v dx \]

Each of these bilinear forms and linear forms lead to matrices and vectors

\[ A_k \qquad\text{and}\qquad \underline f_k. \]

Using the jump matrices \( B_k \), which are provided by the class gsIetiMapper, we obtain the following problem:

\[ \begin{pmatrix} A_1 & & & & B_1^\top \\ & A_2 & & & B_2^\top \\ & &\ddots& & \vdots \\ & & & A_K & B_K^\top \\ B_1 & B_2 &\cdots& B_K & 0 \\ \end{pmatrix} \begin{pmatrix} \underline u_1 \\ \underline u_2 \\ \vdots \\ \underline u_K \\ \underline \lambda \end{pmatrix} = \begin{pmatrix} \underline f_1 \\ \underline f_2 \\ \vdots \\ \underline f_K \\ 0 \end{pmatrix} \]

This is the kind of systems, the class gsIetiSystem is working with. The problem is that – in general – the matrices \( A_k \) are singular. To circumvent this problem, in FETI-dp / IETI-dp methods, we use primal degrees of freedom, like vertex values or edge averages.

These primal dofs are handled by the classes gsPrimalSystem and again gsIetiMapper. The IETI mapper provides vectors \( c_k^{(i)} \) that realize the constraint that the primal dof is set to zero. For each patch, these primal constraints build a matrix \( C_k \). The IETI mapper also provides a vector of indices that associates each of these primal constraints to a (global) index referring to the primal dof.

Finally, we obtain one subspace (the K+1st) that refers to the primal degrees of freedom.

We have the following IETI system:

\[ \begin{pmatrix} \tilde A_1 & & & & \tilde B_1^\top \\ & \tilde A_2 & & & \tilde B_2^\top \\ & & \ddots & & \vdots \\ & & & \tilde A_{K+1} & \tilde B_{K+1}^\top \\ \tilde B_1 & \tilde B_2 & \cdots & \tilde B_{K+1} & 0 \\ \end{pmatrix} \begin{pmatrix} \underline{ \tilde u}_1 \\ \underline{ \tilde u}_2 \\ \vdots \\ \underline{ \tilde u}_{K+1} \\ \underline \lambda \end{pmatrix} = \begin{pmatrix} \underline{ \tilde f}_1 \\ \underline{ \tilde f}_2 \\ \vdots \\ \underline{ \tilde f}_{K+1} \\ 0 \end{pmatrix}, \qquad (*) \]

which is passed to the class gsIetiSystem. For \( k=1,\ldots,K \), we have the local saddle point system

\[ \tilde{A}_k = \begin{pmatrix} A_k & C_k^\top \\ C_k & 0 \end{pmatrix}, \]

the extended jump matrix

\[ \tilde{B}_k = B_k (P_0^{(k)})^\top = \begin{pmatrix} B_k & 0_k \end{pmatrix} \]

and the extended right-hand side

\[ \underline{\tilde f}_k = P_0^{(k)} \underline f_k = \begin{pmatrix} \underline f_k \\ 0 \end{pmatrix}. \]

These modifications are provided by gsPrimalSystem. The class gsPrimalSystem moreover computes the energy minimizing basis functions and represents them as matrices \( \Psi_k = (P_c^{(k)})^\top \tilde{A}_k^{-1} R_c \), where \( R_c=(0,I)^\top \) and \( P_c=(I,0) \). It then computes

\[ \tilde A_{K+1} = \sum_{k=1}^K \Psi_k A_k \Psi_k^\top \quad\text{and}\quad \tilde B_{K+1} = \sum_{k=1}^K B_k \Psi_k^\top \quad\text{and}\quad \underline{\tilde f}_{K+1} = \sum_{k=1}^K \Psi_k \underline f_k. \]

By setting a corresponding flag, it is possible to eliminate the pointwise constraints. A constraint is pointwise the correspoding row of \( C_k \) has only one non-zero entries. Typically, the corner values are pointwise constraints.

Instead of directly solving the IETI-saddle point system (*), we solve the corresponding Schur complement system

\[ F \underline \lambda = \underline g, \]

where

\[ F = \sum_{k=1}^{K+1} \tilde B_k \tilde A_k^{-1} \tilde B_k^\top \quad\text{and}\quad \underline g = \sum_{k=1}^{K+1} \tilde B_k \tilde A_k^{-1} \underline{\tilde f}_k. \]

The Schur complement system is preconditioned with the scaled Dirichlet preconditioner

\[ M_{sD} = \sum_{k=1}^{K} \hat B_k D_k^{-1} S_k D_k^{-1} \hat B_k^\top, \]

where \( S_k \) is the Schur complement of \( A_k \) (not: \(\tilde A_k\)) with respect to the dofs on the interfaces. \( \hat B_k \) is the restriction of \( B_k \) to the dofs on the interfaces. \( D_k \) is a scaling matrix, e.g., based on multiplicity scaling. The class gsScaledDirichletPrec not only collects the contributions of the preconditioner; it also helps with restricting to the skeleton and computing the Schur complement.

In the example file ieti2_example.cpp, we solve the saddle point system (*) with a MINRES solver, preconditioned with a block-diagonal preconditioner, where the last block is \( M_{sD} \).

First, we set up the gsIetiMapper object, which helps us with assembling and keeping track of the dof indices.

    gsIetiMapper<> ietiMapper;

For initializing the IETI mapper, it is necessary to provide a gsDofMapper for the global problem and the function values for the eliminated dofs (typically the Dirichlet conditions). If systems are considered, the method it is expected to have one IETI mapper per variable.

For obtaining a gsDofMapper, we use the assembler:

    {
        typedef gsExprAssembler<>::space  space;
        gsExprAssembler<> assembler;
        space u = assembler.getSpace(mb);
        bc.setGeoMap(mp);
        u.setup(bc, dirichlet::interpolation, 0);
        ietiMapper.init( mb, u.mapper(), u.fixedPart() );
    }

Now we set the primal degrees of freedom. The following commands compute the (vectors that make up the) constraint matrices \( C_k \) and store them in the gsIetiMapper object.

    if (cornersAsPrimals)
        ietiMapper.cornersAsPrimals();
 
    if (edgesAsPrimals)
        ietiMapper.interfaceAveragesAsPrimals(mp,1);
 
    if (facesAsPrimals)
        ietiMapper.interfaceAveragesAsPrimals(mp,2);

The following command computes the jump matrices \( B_k \) and stores them in the gsIetiMapper object.

    ietiMapper.computeJumpMatrices(fullyRedundant, noLagrangeMultipliersForCorners);

Now, we set up the classes for the IETI system, the primal system, and the scaled Dirichlet preconditioner.

    // The ieti system does not have a special treatment for the
    // primal dofs. They are just one more subdomain
    gsIetiSystem<> ieti;
    ieti.reserve(nPatches+1);
 
    // The scaled Dirichlet preconditioner is independent of the
    // primal dofs.
    gsScaledDirichletPrec<> prec;
    prec.reserve(nPatches);
 
    // Setup the primal system, which needs to know the number of primal dofs.
    gsPrimalSystem<> primal(ietiMapper.nPrimalDofs());
    if (eliminateCorners)
        primal.setEliminatePointwiseConstraints(true);

Now, we assemble the linear system and the right-hand side on each patch using a for loop. We use local versions of the geometry (gsMultiPatch), the basis (gsMultiBasis) and the boundary conditions (gsBoundaryConditions).

The function initFeSpace provides a new dof mapper and the Dirichlet data. This is necessary since it might happen that a 2d-patch touches the Dirichlet boundary just with a corner or that a 3d-patch touches the Dirichlet boundary with a corner or an edge. These cases are not covered by bc.getConditionsForPatch

The example file ieti2_example.cpp shows how to assemble with a classical assembler.

    for (index_t k=0; k<nPatches; ++k)
    {
        // We use the local variants of everything
        gsBoundaryConditions<> bc_local;
        bc.getConditionsForPatch(k,bc_local);
        gsMultiPatch<> mp_local = mp[k];
        gsMultiBasis<> mb_local = mb[k];
 
        // The usual stuff for the expression assembler
        typedef gsExprAssembler<>::geometryMap geometryMap;
        typedef gsExprAssembler<>::variable    variable;
        typedef gsExprAssembler<>::space       space;
        typedef gsExprAssembler<>::solution    solution;
 
        // We set up the assembler
        gsExprAssembler<> assembler(1,1);
 
        // Elements used for numerical integration
        assembler.setIntegrationElements(mb_local);
        gsExprEvaluator<> ev(assembler);
 
        // Set the geometry map
        geometryMap G = assembler.getMap(mp_local);
 
        // Set the discretization space
        space u = assembler.getSpace(mb_local);
 
        // Incorporate Dirichlet BC
        bc_local.setGeoMap(mp_local);
        u.setup(bc_local, dirichlet::interpolation, 0);
 
        // This function provides a new dof mapper and the Dirichlet data
        // This is necessary since it might happen that a 2d-patch touches the
        // Dirichlet boundary just with a corner or that a 3d-patch touches the
        // Dirichlet boundary with a corner or an edge. These cases are not
        // covered by bc.getConditionsForPatch
        ietiMapper.initFeSpace(u,k);
 
        // Set the source term
        auto ff = assembler.getCoeff(f, G);
 
        // Initialize the system
        assembler.initSystem();
 
        // Compute the system matrix and right-hand side
        assembler.assemble( igrad(u, G) * igrad(u, G).tr() * meas(G), u * ff * meas(G) );
 
        // Add contributions from Neumann conditions to right-hand side
        variable g_N = assembler.getBdrFunction();
        assembler.assembleBdr(bc_local.get("Neumann"),  u * g_N.val() * nv(G).norm() );
 
        // Fetch data
        gsSparseMatrix<real_t, RowMajor> jumpMatrix  = ietiMapper.jumpMatrix(k);
        gsSparseMatrix<>                 localMatrix = assembler.matrix();
        gsMatrix<>                       localRhs    = assembler.rhs();

Now, we tell the preconditioner about the local matrix \( A_k \), the local vector \( \underline f_k \) and the matrix \( B_k \). This is done before the class gsPrimalSystem modifies these objects. The following command does all necessary steps:

        prec.addSubdomain(
            gsScaledDirichletPrec<>::restrictToSkeleton(
                jumpMatrix,
                localMatrix,
                ietiMapper.skeletonDofs(k)
            )
        );

This can also be done in a step-by-step way, for example if one does not want to use a sparse Cholesky solver (which requires symmetric positive definite problems and would be done by default). This is shown in ieti2_example.cpp.

As next step, we treat the primal dofs. We first compute energy minimizing primal basis functions. Then the following commands amends the variables jumpMatrix, localMatrix, and localRhs, which hold \( B_k \), \( A_k \) and \( \underline f_k \) before calling the function and \( \tilde B_k \), \( \tilde A_k \), and \( \underline{\tilde f}_k \) thereafter. Finally, this function collects the contributions for \( \tilde B_{K+1} \), \( \tilde A_{K+1} \), and \( \underline{\tilde f}_{K+1} \) in the object.

        primal.handleConstraints(
            ietiMapper.primalConstraints(k),
            ietiMapper.primalDofIndices(k),
            jumpMatrix,
            localMatrix,
            localRhs
        );

This can also be done in a step-by-step way, for example if one does not want to use a sparse LU solver (which would be done by default). This is shown in the file ieti2_example.cpp.

Now, we add \( \tilde B_k \), \( \tilde A_k \), and \( \underline{\tilde f}_k \) to the IETI-system:

        ieti.addSubdomain(
            jumpMatrix.moveToPtr(),
            makeMatrixOp(localMatrix.moveToPtr()),
            give(localRhs)
        );

Additionally, one can also provide a solver object that realizes \( \tilde A_k^{-1} \), possibly the one which has already been used for the setup of the bases for the primal doefs. By default, again, a LU solver would be used.

    } // end for

Finally, we add \( \tilde B_{K+1} \), \( \tilde A_{K+1} \), and \( \underline{\tilde f}_{K+1} \) to the IETI-system:

    if (ietiMapper.nPrimalDofs()>0)
    {
        // It is not required to provide a local solver to .addSubdomain,
        // since a sparse LU solver would be set up on the fly if required.
        // Here, we make use of the fact that we can use a Cholesky solver
        // because the primal problem is symmetric and positive definite:
        gsLinearOperator<>::Ptr localSolver
            = makeSparseCholeskySolver(primal.localMatrix());
 
        ieti.addSubdomain(
            primal.jumpMatrix().moveToPtr(),
            makeMatrixOp(primal.localMatrix().moveToPtr()),
            give(primal.localRhs()),
            localSolver
        );
    }

The gsScaledDirichletPrec can automatically compute the multiplicity scaling. Other scaling matrices can be manually defined.

    prec.setupMultiplicityScaling();

Now, we compute the right-hand-side

\[ \underline g = \sum_{k=1}^{K+1} \tilde B_k \tilde A_k^{-1} \underline{\tilde f}_k. \]

    gsMatrix<> rhsForSchur = ieti.rhsForSchurComplement();

We choose a random initial guess.

    gsMatrix<> lambda;
    lambda.setRandom( ieti.nLagrangeMultipliers(), 1 );

As a next step, we solve the problem

\[ F \underline \lambda = \underline g \]

using the scaled Dirichlet preconditioner with a conjugate gradient solver.

    gsConjugateGradient<> PCG( ieti.schurComplement(), prec.preconditioner() );
    PCG.setOptions( cmd.getGroup("Solver") ).solveDetailed( rhsForSchur, lambda, errorHistory );

So far, we have only computed the Lagrange multipliers \( \underline \lambda \), from which we want to recover the solution.

Here, we first compute the solutions \( \underline{\tilde u}_k \) for \( k=1,\ldots,K+1 \). Here, the first K solutions refer to the patches and have homogenous values for the primal dofs. So, we distribute the primal dof values from the primal problem (= last subdomain) to the patches (=first K subdomains) and obtain the solutions \( \underline u_k \) for \( k=1,\ldots,K \). Then, finally, the IETI mapper is able to combine everything into one solution vector \( \underline u \).

    std::vector<gsMatrix<>> uLocal = primal.distributePrimalSolution(
        ieti.constructSolutionFromLagrangeMultipliers(lambda)
    );
    gsMatrix<> uGlobal = ietiMapper.constructGlobalSolutionFromLocalSolutions(uLocal);

Annotated source file

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

 
#include <ctime>
#include <gismo.h>
 
using namespace gismo;
 
int main(int argc, char *argv[])
{
    /************** Define command line options *************/
 
    std::string geometry("domain2d/yeti_mp2.xml");
    index_t splitPatches = 1;
    real_t stretchGeometry = 1;
    index_t refinements = 1;
    index_t degree = 2;
    std::string boundaryConditions("d");
    std::string primals("c");
    bool eliminateCorners = false;
    real_t tolerance = 1.e-8;
    index_t maxIterations = 100;
    bool calcEigenvalues = false;
    std::string out;
    bool plot = false;
 
    gsCmdLine cmd("Solves a PDE with an isogeometric discretization using an isogeometric tearing and interconnecting (IETI) solver.");
    cmd.addString("g", "Geometry",              "Geometry file", geometry);
    cmd.addInt   ("",  "SplitPatches",          "Split every patch that many times in 2^d patches", splitPatches);
    cmd.addReal  ("",  "StretchGeometry",       "Stretch geometry in x-direction by the given factor", stretchGeometry);
    cmd.addInt   ("r", "Refinements",           "Number of uniform h-refinement steps to perform before solving", refinements);
    cmd.addInt   ("p", "Degree",                "Degree of the B-spline discretization space", degree);
    cmd.addString("b", "BoundaryConditions",    "Boundary conditions", boundaryConditions);
    cmd.addString("c", "Primals",               "Primal constraints (c=corners, e=edges, f=faces)", primals);
    cmd.addSwitch("e", "EliminateCorners",      "Eliminate corners (if they are primals)", eliminateCorners);
    cmd.addReal  ("t", "Solver.Tolerance",      "Stopping criterion for linear solver", tolerance);
    cmd.addInt   ("",  "Solver.MaxIterations",  "Maximum iterations for linear solver", maxIterations);
    cmd.addSwitch("",  "Solver.CalcEigenvalues","Estimate eigenvalues based on Lanczos", calcEigenvalues);
    cmd.addString("",  "out",                   "Write solution and used options to file", out);
    cmd.addSwitch(     "plot",                  "Plot the result with Paraview", plot);
 
    try { cmd.getValues(argc,argv); } catch (int rv) { return rv; }
 
 
    if ( ! gsFileManager::fileExists(geometry) )
    {
        gsInfo << "Geometry file could not be found.\n";
        gsInfo << "I was searching in the current directory and in: " << gsFileManager::getSearchPaths() << "\n";
        return EXIT_FAILURE;
    }
 
    gsInfo << "Run ieti_example with options:\n" << cmd << std::endl;
 
    /******************* Define geometry ********************/
 
    gsInfo << "Define geometry... " << std::flush;
 
    gsMultiPatch<>::uPtr mpPtr = gsReadFile<>(geometry);
    if (!mpPtr)
    {
        gsInfo << "No geometry found in file " << geometry << ".\n";
        return EXIT_FAILURE;
    }
    gsMultiPatch<>& mp = *mpPtr;
 
    for (index_t i=0; i<splitPatches; ++i)
    {
        gsInfo << "split patches uniformly... " << std::flush;
        mp = mp.uniformSplit();
    }
 
    if (stretchGeometry!=1)
    {
       gsInfo << "and stretch it... " << std::flush;
       for (size_t i=0; i!=mp.nPatches(); ++i)
           const_cast<gsGeometry<>&>(mp[i]).scale(stretchGeometry,0);
       // Const cast is allowed since the object itself is not const. Stretching the
       // overall domain keeps its topology.
    }
 
    gsInfo << "done.\n";
 
    /************** Define boundary conditions **************/
 
    gsInfo << "Define right-hand-side and boundary conditions... " << std::flush;
 
    // Right-hand-side
    gsFunctionExpr<> f( "2*sin(x)*cos(y)", mp.geoDim() );
 
    // Dirichlet function
    gsFunctionExpr<> gD( "sin(x)*cos(y)", mp.geoDim() );
 
    // Neumann
    gsConstantFunction<> gN( 1.0, mp.geoDim() );
 
    gsBoundaryConditions<> bc;
    {
        const index_t len = boundaryConditions.length();
        index_t i = 0;
        for (gsMultiPatch<>::const_biterator it = mp.bBegin(); it < mp.bEnd(); ++it)
        {
            char b_local;
            if ( len == 1 )
                b_local = boundaryConditions[0];
            else if ( i < len )
                b_local = boundaryConditions[i];
            else
            {
                gsInfo << "\nNot enough boundary conditions given.\n";
                return EXIT_FAILURE;
            }
 
            if ( b_local == 'd' )
                bc.addCondition( *it, condition_type::dirichlet, &gD );
            else if ( b_local == 'n' )
                bc.addCondition( *it, condition_type::neumann, &gN );
            else
            {
                gsInfo << "\nInvalid boundary condition given; only 'd' (Dirichlet) and 'n' (Neumann) are supported.\n";
                return EXIT_FAILURE;
            }
 
            ++i;
        }
        if ( len > i )
            gsInfo << "\nToo many boundary conditions have been specified. Ignoring the remaining ones.\n";
        gsInfo << "done. "<<i<<" boundary conditions set.\n";
    }
 
 
    /************ Setup bases and adjust degree *************/
 
    gsMultiBasis<> mb(mp);
 
    gsInfo << "Setup bases and adjust degree... " << std::flush;
 
    for ( size_t i = 0; i < mb.nBases(); ++ i )
        mb[i].setDegreePreservingMultiplicity(degree);
 
    for ( index_t i = 0; i < refinements; ++i )
        mb.uniformRefine();
 
    gsInfo << "done.\n";
 
    /********* Setup assembler and assemble matrix **********/
 
    gsInfo << "Setup assembler and assemble matrix... " << std::flush;
 
    const index_t nPatches = mp.nPatches();
 
    gsIetiMapper<> ietiMapper;
 
    // We start by setting up a global FeSpace that allows us to
    // obtain a dof mapper and the Dirichlet data
    {
        typedef gsExprAssembler<>::space  space;
        gsExprAssembler<> assembler;
        space u = assembler.getSpace(mb);
        bc.setGeoMap(mp);
        u.setup(bc, dirichlet::interpolation, 0);
        ietiMapper.init( mb, u.mapper(), u.fixedPart() );
    }
 
    // Which primal dofs should we choose?
    bool cornersAsPrimals = false, edgesAsPrimals = false, facesAsPrimals = false;
    for (size_t i=0; i<primals.length(); ++i)
        switch (primals[i])
        {
            case 'c': cornersAsPrimals = true;   break;
            case 'e': edgesAsPrimals = true;     break;
            case 'f': facesAsPrimals = true;     break;
            default:
                gsInfo << "\nUnkown type of primal constraint: \"" << primals[i] << "\"\n";
                return EXIT_FAILURE;
        }
 
    // We tell the ieti mapper which primal constraints we want; calling
    // more than one such function is possible.
    if (cornersAsPrimals)
        ietiMapper.cornersAsPrimals();
 
    if (edgesAsPrimals)
        ietiMapper.interfaceAveragesAsPrimals(mp,1);
 
    if (facesAsPrimals)
        ietiMapper.interfaceAveragesAsPrimals(mp,2);
 
    // Compute the jump matrices
    bool fullyRedundant = true,
         noLagrangeMultipliersForCorners = cornersAsPrimals;
    ietiMapper.computeJumpMatrices(fullyRedundant, noLagrangeMultipliersForCorners);
 
    // The ieti system does not have a special treatment for the
    // primal dofs. They are just one more subdomain
    gsIetiSystem<> ieti;
    ieti.reserve(nPatches+1);
 
    // The scaled Dirichlet preconditioner is independent of the
    // primal dofs.
    gsScaledDirichletPrec<> prec;
    prec.reserve(nPatches);
 
    // Setup the primal system, which needs to know the number of primal dofs.
    gsPrimalSystem<> primal(ietiMapper.nPrimalDofs());
    if (eliminateCorners)
        primal.setEliminatePointwiseConstraints(true);
 
    for (index_t k=0; k<nPatches; ++k)
    {
        // We use the local variants of everything
        gsBoundaryConditions<> bc_local;
        bc.getConditionsForPatch(k,bc_local);
        gsMultiPatch<> mp_local = mp[k];
        gsMultiBasis<> mb_local = mb[k];
 
        // The usual stuff for the expression assembler
        typedef gsExprAssembler<>::geometryMap geometryMap;
        typedef gsExprAssembler<>::variable    variable;
        typedef gsExprAssembler<>::space       space;
        typedef gsExprAssembler<>::solution    solution;
 
        // We set up the assembler
        gsExprAssembler<> assembler(1,1);
 
        // Elements used for numerical integration
        assembler.setIntegrationElements(mb_local);
        gsExprEvaluator<> ev(assembler);
 
        // Set the geometry map
        geometryMap G = assembler.getMap(mp_local);
 
        // Set the discretization space
        space u = assembler.getSpace(mb_local);
 
        // Incorporate Dirichlet BC
        bc_local.setGeoMap(mp_local);
        u.setup(bc_local, dirichlet::interpolation, 0);
 
        // This function provides a new dof mapper and the Dirichlet data
        // This is necessary since it might happen that a 2d-patch touches the
        // Dirichlet boundary just with a corner or that a 3d-patch touches the
        // Dirichlet boundary with a corner or an edge. These cases are not
        // covered by bc.getConditionsForPatch
        ietiMapper.initFeSpace(u,k);
 
        // Set the source term
        auto ff = assembler.getCoeff(f, G);
 
        // Initialize the system
        assembler.initSystem();
 
        // Compute the system matrix and right-hand side
        assembler.assemble( igrad(u, G) * igrad(u, G).tr() * meas(G), u * ff * meas(G) );
 
        // Add contributions from Neumann conditions to right-hand side
        variable g_N = assembler.getBdrFunction();
        assembler.assembleBdr(bc_local.get("Neumann"),  u * g_N.val() * nv(G).norm() );
 
        // Fetch data
        gsSparseMatrix<real_t, RowMajor> jumpMatrix  = ietiMapper.jumpMatrix(k);
        gsSparseMatrix<>                 localMatrix = assembler.matrix();
        gsMatrix<>                       localRhs    = assembler.rhs();
 
        // Add the patch to the scaled Dirichlet preconditioner
        prec.addSubdomain(
            gsScaledDirichletPrec<>::restrictToSkeleton(
                jumpMatrix,
                localMatrix,
                ietiMapper.skeletonDofs(k)
            )
        );
 
        // This function writes back to jumpMatrix, localMatrix, and localRhs,
        // so it must be called after prec.addSubdomain().
        primal.handleConstraints(
            ietiMapper.primalConstraints(k),
            ietiMapper.primalDofIndices(k),
            jumpMatrix,
            localMatrix,
            localRhs
        );
 
        // Add the patch to the Ieti system
        ieti.addSubdomain(
            jumpMatrix.moveToPtr(),
            makeMatrixOp(localMatrix.moveToPtr()),
            give(localRhs)
        );
    } // end for
 
    // Add the primal problem if there are primal constraints
    if (ietiMapper.nPrimalDofs()>0)
    {
        // It is not required to provide a local solver to .addSubdomain,
        // since a sparse LU solver would be set up on the fly if required.
        // Here, we make use of the fact that we can use a Cholesky solver
        // because the primal problem is symmetric and positive definite:
        gsLinearOperator<>::Ptr localSolver
            = makeSparseCholeskySolver(primal.localMatrix());
 
        ieti.addSubdomain(
            primal.jumpMatrix().moveToPtr(),
            makeMatrixOp(primal.localMatrix().moveToPtr()),
            give(primal.localRhs()),
            localSolver
        );
    }
 
    gsInfo << "done. " << ietiMapper.nPrimalDofs() << " primal dofs.\n";
 
    /**************** Setup solver and solve ****************/
 
    gsInfo << "Setup solver and solve... \n"
        "    Setup multiplicity scaling... " << std::flush;
 
    // Tell the preconditioner to set up the scaling
    prec.setupMultiplicityScaling();
 
    gsInfo << "done.\n    Setup rhs... " << std::flush;
    // Compute the Schur-complement contribution for the right-hand-side
    gsMatrix<> rhsForSchur = ieti.rhsForSchurComplement();
 
    gsInfo << "done.\n    Setup cg solver for Lagrange multipliers and solve... " << std::flush;
    // Initial guess
    gsMatrix<> lambda;
    lambda.setRandom( ieti.nLagrangeMultipliers(), 1 );
 
    gsMatrix<> errorHistory;
 
    // This is the main cg iteration
    gsConjugateGradient<> PCG( ieti.schurComplement(), prec.preconditioner() );
    PCG.setOptions( cmd.getGroup("Solver") ).solveDetailed( rhsForSchur, lambda, errorHistory );
 
    gsInfo << "done.\n    Reconstruct solution from Lagrange multipliers... " << std::flush;
    // Now, we want to have the global solution for u
    std::vector<gsMatrix<>> uLocal = primal.distributePrimalSolution(
        ieti.constructSolutionFromLagrangeMultipliers(lambda)
    );
    gsMatrix<> uGlobal = ietiMapper.constructGlobalSolutionFromLocalSolutions(uLocal);
    gsInfo << "done.\n\n";
 
    /******************** Print end Exit ********************/
 
    const index_t iter = errorHistory.rows()-1;
    const bool success = errorHistory(iter,0) < tolerance;
    if (success)
        gsInfo << "Reached desired tolerance after " << iter << " iterations:\n";
    else
        gsInfo << "Did not reach desired tolerance after " << iter << " iterations:\n";
 
    if (errorHistory.rows() < 20)
        gsInfo << errorHistory.transpose() << "\n\n";
    else
        gsInfo << errorHistory.topRows(5).transpose() << " ... " << errorHistory.bottomRows(5).transpose()  << "\n\n";
 
    if (calcEigenvalues)
        gsInfo << "Estimated condition number: " << PCG.getConditionNumber() << "\n";
 
    if (!out.empty())
    {
        gsFileData<> fd;
        std::time_t time = std::time(NULL);
        fd.add(cmd);
        fd.add(uGlobal);
        fd.addComment(std::string("ieti_example   Timestamp:")+std::ctime(&time));
        fd.save(out);
        gsInfo << "Write solution to file " << out << "\n";
    }
 
    if (plot)
    {
        gsInfo << "Write Paraview data to file ieti_result.pvd\n";
        gsMultiPatch<> mpsol;
        for (index_t k=0; k<nPatches; ++k)
            mpsol.addPatch( mb[k].makeGeometry( ietiMapper.incorporateFixedPart(k, uLocal[k])  ) );
        gsWriteParaview<>( gsField<>( mp, mpsol ), "ieti_result", 1000);
        //gsFileManager::open("ieti_result.pvd");
    }
 
    if (!plot&&out.empty())
    {
        gsInfo << "Done. No output created, re-run with --plot to get a ParaView "
                  "file containing the solution or --out to write solution to xml file.\n";
    }
    return success ? EXIT_SUCCESS : EXIT_FAILURE;
}