Back scattering from a PEC disk with lowest order BGT boundary condition

---

Defining the mesh

Option 1: Reading mesh from files

The mesh will be imported from four files in the usual format in Concepts.

       concepts::Import2dMeshGeneral msh
       (SOURCEDIR "/app-mengyu/PEC-diskCoord.dat",
       SOURCEDIR "/app-mengyu/PEC-diskElements.dat",
       concepts::AttributesFile(SOURCEDIR "/app-mengyu/PEC-diskBoundary.dat"),
       concepts::EdgRadiaFile(SOURCEDIR "/app-mengyu/PEC-diskEdgRadia.dat"));

Option 2: Generating mesh through class concepts::Circle

Alternatively, for the circular geometry, the mesh can also be imported using class concepts::Circle.

       Real innerRadius(1.0);
       Real R(8.0); // the outer radius 
       Real rdata[4] = {1.5,2.0,4.0,R};
       concepts::Array<Real> rings(rdata,4);
       concepts::Circle msh(innerRadius,rings);

---

The attributes of the boundary edges will be assigned during mesh generation.

---

The mesh data can be written to a MATLAB file, in order to view the mesh. About the output format, please refer to Output of FE solution on 2D hp-adaptive space .

       graphics::MatlabGraphics(msh, "PEC_disk.m");

Build the Space

The space will be built with class hp2D::hpAdaptiveSpaceH1. The trace space will be built with class hp2D::TraceSpace.
The inner trace space tspcInner will be built according to its attribute 2, and the outer tspcOuter will be built according to its attribute 1.

       uint p = 3;  // polynomial degree
       uint l = 1;   // number of refinement
       hp2D::hpAdaptiveSpaceH1 space(msh, l, p);
       space.rebuild();
       hp2D::TraceSpace tspcInner(space, concepts::makeSet<uint>(1, 2),
                               hp2D::TraceSpace::FIRST);
       hp2D::TraceSpace tspcOuter(space, concepts::makeSet<uint>(1, 1),
                               hp2D::TraceSpace::FIRST);

Construct right hand side

Parse the formula

The complex formula will be parsed using the constructor of class concepts::ParsedFormula.

    sR << "((" << k << ")"<<"*x/sqrt(x*x+y*y)*sin(x*"<<"("<<k<<")"<<"))";
    sI << "((-1)*(" << k << ")"<<"*x/sqrt(x*x+y*y)*cos(x*"<<"("<<k<<")"<<"))"; 
    concepts::ParsedFormula<Cmplx> frm(sR.str(), sI.str());

Construct linear form

linear form lform is used to construct right hand side, to obtain lform on trace space, class hp1D::Riesz is applied.

    hp1D::Riesz<Cmplx> lform(frm);

Construct right hand side

Finally, the right hand side will be constructed by the constructor of class concepts::Vector.

    concepts::Vector<Cmplx> rhs(tspcInner, lform);

Build matrices

Matrices building

The matrices are built by the constructor of class concepts::SparseMatrix.

    concepts::SparseMatrix<Cmplx> S(space.dim(), space.dim());

will give the empty sparse matrix, and

    hp2D::Laplace<> la;
    concepts::SparseMatrix<Real> stiff(space, la);

will construct matrices by bilinear form.

Summing up the matrices

for sparse matrices A and B, the function A.add(B,coeff) adds A itself to B with multiplier coeff. i.e. B <- B + coeff*A

    stiff.addInto(S, 1.0);    // stiff adds itself to S

Solve the system

The linear system is solved with a sparse direct solver SuperLU

    concepts::SuperLU<Cmplx> solver(S);
    concepts::Vector<Cmplx> sol(space);   // coefficient vector of the solution
    solver(rhs, sol);

Output solution

The solution will be exported to a MATLAB file. To avoid losing precision with interpolation, recomputeShapefunctions() is developed to compute precisely.

    hp2D::IntegrableQuad::rule().set(concepts::TRAPEZE, true, 12);
    space.recomputeShapefunctions();        // recomputation on new quadrature points
    graphics::MatlabGraphics(space, "PEC_wedge_sol.m", sol);