strong form finite element formulation prototype

ex0strong.cpp

//                                MFEM Example 0
//
// Compile with: make ex0
//
// Sample runs:  ex0
//               ex0 -m ../data/fichera.mesh
//               ex0 -m ../data/square-disc.mesh -o 2
//
// Description: This example code demonstrates the most basic usage of MFEM to
//              define a simple finite element discretization of the Laplace
//              problem -Delta u = 1 with zero Dirichlet boundary conditions.
//              General 2D/3D mesh files and finite element polynomial degrees
//              can be specified by command line options.

#include "mfem.hpp"
#include <fstream>
#include <iostream>

using namespace std;
using namespace mfem;

class StrongDiffusionIntegrator : public DiffusionIntegrator {
public:
    int flag = 0;
    void AssembleElementMatrix(const FiniteElement &el,
                               ElementTransformation &Trans,
                               DenseMatrix &elmat) override {

        int dim = el.GetDim();
        int ndof = el.GetDof();
        const IntegrationRule *ir = IntRule ? IntRule : &GetRule(el, el);
        int n_gauss = ir->GetNPoints();
        Vector shape(ndof); // shape function
        Vector d2shape(ndof); // shape function second derivative

        DenseMatrix dshape(ndof, dim);
        DenseMatrix dshape_eps(ndof, dim);
        elmat.SetSize(ndof);
        elmat = 0.0;
        const real_t eps = 1e-6;

        if (!flag) {
            int n = 100;
            double step = 1.0 / (n - 1);
            ofstream ofs("phi.dat");
            for (int i = 0; i < n; i++) {
                IntegrationPoint ip = {i * step, 0, 0, 0};
                el.CalcShape(ip, shape);
                ofs << ip.x << " ";
                for (int j = 0; j < ndof; j++) {
                    ofs << shape(j) << " ";
                }
                ofs << endl;
            }
        }

        for (int q = 0; q < ir->GetNPoints(); q++) {
            const IntegrationPoint &ip = ir->IntPoint(q);
            real_t w = ip.weight;
            el.CalcShape(ip, shape);

            // use finite difference to calculate the second derivative
            // assume 1D
            IntegrationPoint ip2 = ip;
            ip2.x += eps;
            el.CalcDShape(ip, dshape);
            el.CalcDShape(ip2, dshape_eps);
            for (int i = 0; i < ndof; i++) {
                d2shape(i) = (dshape_eps(i, 0) - dshape(i, 0)) / eps;
            }

            if (!flag) {
                for (int i = 0; i < ndof; i++) {
                    std::printf("IP %10f %10f %10f %10f\n", ip.x, shape(i), dshape(i, 0), d2shape(i));
                }
            }

            for (int i = 0; i < ndof; i++) {
                for (int j = 0; j < ndof; j++) {
                    // Weak form (solution is correct)
                    //elmat(i, j) += dshape(i, 0) * dshape(j, 0) *  w * Trans.Weight();

                    // Strong form (solution is wrong)
                    elmat(i, j) += d2shape(j) * shape(i) * w * Trans.Weight();
                }
            }
        }
        if (!flag) {
            elmat.Print();
            DenseMatrix tmp(ndof);
            tmp = elmat;
            tmp.Transpose();
            tmp.Add(-1.0, elmat);
            real_t max_val = tmp.MaxMaxNorm();
            std::cout << "||K - K^t|| " << max_val << std::endl;
            if (max_val > 1e-10) {
                std::cout << "ERROR: K is not symmetric" << std::endl;
            }
        }

        flag = 1;
    }
};

int main(int argc, char *argv[]) {
    int iterative = 0;
    int order = 2;
    Mesh mesh(Mesh::MakeCartesian1D(5, 1));

    // 3. Define a finite element space on the mesh. Here we use H1 continuous
    //    high-order Lagrange finite elements of the given order.
    H1_FECollection fec(order, mesh.Dimension());
    FiniteElementSpace fespace(&mesh, &fec);
    cout << "Number of unknowns: " << fespace.GetTrueVSize() << endl;

    // 4. Extract the list of all the boundary DOFs. These will be marked as
    //    Dirichlet in order to enforce zero boundary conditions.
    Array<int> boundary_dofs;
    fespace.GetBoundaryTrueDofs(boundary_dofs);

    // 5. Define the solution x as a finite element grid function in fespace. Set
    //    the initial guess to zero, which also sets the boundary conditions.
    GridFunction x(&fespace);
    x = 0.0;

    // 6. Set up the linear form b(.) corresponding to the right-hand side.
    ConstantCoefficient one(1.0);
    LinearForm b(&fespace);
    b.AddDomainIntegrator(new DomainLFIntegrator(one));
    b.Assemble();

    // 7. Set up the bilinear form a(.,.) corresponding to the -Delta operator.
    BilinearForm a(&fespace);
    a.AddDomainIntegrator(new StrongDiffusionIntegrator);
    a.Assemble();

    // 8. Form the linear system A X = B. This includes eliminating boundary
    //    conditions, applying AMR constraints, and other transformations.
    SparseMatrix A;
    Vector B, X;
    a.FormLinearSystem(boundary_dofs, x, b, A, X, B);

    // 9. Solve the system
    // NOTICE: because the system is not symmetric, the PCG solver fails
    if (iterative) {
        GSSmoother M(A);
        PCG(A, M, B, X, 1, 200, 1e-12, 0.0);
    } else {
        UMFPackSolver umf_solver;
        umf_solver.SetOperator(A);
        umf_solver.Mult(B, X);
    }

    // 10. Recover the solution x as a grid function and save to file. The output
    //     can be viewed using GLVis as follows: "glvis -m mesh.mesh -g sol.gf"
    a.RecoverFEMSolution(X, b, x);

    x.Save("sol.gf");
    mesh.Save("mesh.mesh");

    return 0;
}