Overview

Finite element and virtual element library with support for polytopic elements.

Building

Requirements

In order to build PTEMS, the following is required:

• CMake 3.13 (or later)
• CMake-compatible build system such as Make, Ninja, or Visual Studio^[MSVC]
• C++17 compiler, depending on platform and build system:
  • g++ 8.0.0 or later
  • clang++ 8.0.0^[AppleClang] or later (TODO: May need 11 if use shared_ptr to arrays)
  • Visual Studio 2019^[MSVC] or later (either the standard MSVC compiler, or the clang compiler)

PTEMS is currently supported on Linux, expected to work on macOS but untested, and partially supported (without MUMPS) on Windows under 32-bit builds ONLY, cf. issue #2. Windows support on Windows Subsystem for Linux, MSYS, or similar, is untested.

Basic Build Instructions

The following instructions document how to build PTEMS' with default options using the current configured or default CMake generator and C++ compiler for the platform. For details on how to set the generator/compiler to use see the CMake documentation.

We strongly suggest reading the rest of the build instructions before using the commands in this section.

Note, these instructions are based on the use of the standard command line in Linux/macOS or a suitable developer command prompt on Windows (such as the Visual Studio Developer Command Prompt), and require that the cmake command is available on the search PATH. Alternatively, the project folder can be used from Visual Studio Code, providing C++ and CMake Tools extensions are installed; then, you can configure and build using Visual Studio Code's command palette, see Get started with CMake Tools.

Change to the directory containing the PTEMS source and then configure PTEMS with the following command:

cmake -B build

See Build Options below for more detailed options to configure PTEMS.

Then build PTEMS by executing the following command:

cmake --build build

Note that the build system can perform multiple 'jobs' at once (by running multiple compile processes). CMake and the build system will choice a default for how many jobs to run at once, which may be one or potentially the number of processors on the machine. However, you can specify the number of simultaneous jobs to run with the -j argument

cmake --build build -j N

where N is the number of jobs to run.

Note that running the build command multiple times only compiles changed files (should be none!). If you want to force a rebuild of PTEMS library you can specify --clean-first:

cmake --build build --clean-first

This is NOT necessary the first time build is ran.

You can just clean the build using the clean target:

cmake --build build --target clean

Optional Libraries

The PTEMS library has support for optional external libraries. The basic build instructions above will attempt to do the sensible choice depending on the presence of these libraries (see details for each library below for what the sensible thing is).

MPI

PTEMS contains experimental support for parallel processing using any available C MPI library. PTEMS will by default be built without MPI, but it can be forced to be enabled (see Build Options).

MUMPS

The default linear solver is to use MUMPS (MUltifrontal Massively Parallel sparse direct Solver). If MUMPS is detected by the CMake script then PTMES will be built with MUMPS support, otherwise the default libnear solver will fallback to a built-in GMRES with ILU0 preconditioner (optimiality of solver not guaranteed).

MUMPS is usually as multiple binaries, which each binary being a combination of s data type of the matrix (float, double, complex float, and complex double) and whether the solver is a parallel or sequential solver. Multiple binaries of different data types can be included in PTEMS, but it is not possible to combine parallel and sequential solvers. As such, PTEMS will bind use the parallel versions if MPI is enabled and the sequential version otherwise.

Matrix Data Type	Sequential Library	Parallel Library	Header File	Build Option Name
float	smumps_seq	smumps	smumps_c.h	PTEMS_USE_MUMPS
double	dmumps_seq	dmumps	dmumps_c.h	PTEMS_USE_MUMPS_F
complex float	cmumps_seq	cmumps	cmumps_c.h	PTEMS_USE_MUMPS_COMPLEX
complex double	zmumps_seq	zmumps	zmumps_c.h	PTEMS_USE_MUMPS_COMPLEX_F

The double and complex double variants installed in standard search locations should be found by default. The float and complex float variants are by default not included. Setting the specified build option in the above table to ON, OFF, or DEFAULT, (see Build Options) enables, disables, or enables only if the library is found, respectively, for the matching data type. See Build Options for how to aid the detection.

MUMPS can be built from source but is also available for some Linux distributions via their package manager:

Debian/Ubuntu : lib-mumps-seq-dev and lib-mumps-dev for the sequential and parallel versions for all data types, respectively,

Red Hat/Fedora : MUMPS-devel

METIS

The default mesh agglomeration technique uses METIS. If METIS is detected by the CMake script then PTEMS will be built with METIS support, otherwise without. This opposite option can be forced (see Build Options).

METIS can be built from source but is also available for some Linux distributions via their package manager:

Debian/Ubuntu : libmetis-dev

Red Hat/Fedora : metis-devel

A library installed in standard search locations should be found by default. Similarly if the CMake projects for METIS is added before the project for PTEMS in a parent CMake build it should be detected. See Build Options for how to aid the detection.

Build Options

Several options can be passed to the initial configuration step using the -D command line option

cmake -B build -D NAME=VALUE

where NAME is the name of the option in the below table and VALUE is one of the values listed for the option. If not specified the default value in the below table is used.

Option Name	Values	Default	Description
PTEMS_USE_MPI	ON, OFF	OFF	Build with MPI support
PTEMS_USE_TRIANGLE	ON, OFF	ON	Build with Triangle library, see Included Libraries
PTEMS_TRIANGLE_INCDIR	Folder path	""	Path to directory containing triangle.h from Triangle, see Included Libraries
PTEMS_USE_TETGEN	ON, OFF	ON	Build with TetGen library, see Included Libraries
PTEMS_TETGEN_INCDIR	Folder path	""	Path to directory containing tetgen.h from TetGen, see Included Libraries
PTEMS_USE_METIS	DEFAULT, ON, OFF	DEFAULT	Build with METIS support for mesh agglomeration
PTEMS_USE_MUMPS	DEFAULT, ON, OFF	DEFAULT	Build with support for MUMPS with double precision arithmetic
PTEMS_USE_MUMPS_F	DEFAULT, ON, OFF	OFF	Build with support for MUMPS with single precision arithmetic
PTEMS_USE_MUMPS_COMPLEX	DEFAULT, ON, OFF	DEFAULT	Build with support for MUMPS with double precision arithmetic for complex number
PTEMS_USE_MUMPS_COMPLEX_F	DEFAULT, ON, OFF	OFF	Build with support for MUMPS with single precision arithmetic for complex number
PTEMS_BUILD_TESTS	ON, OFF	ON	Build tests
PTEMS_BUILD_EXAMPLES	ON, OFF	ON	Build the example programs in example
PTEMS_LIBRARY_PATHS	List [List]	""	Specify extra search directories for finding libraries
PTEMS_INCLUDE_PATHS	List [List]	""	Specify extra include directories for header files

ON and OFF enable, or disable, the relevant option. DEFAULT indicates it assumes ON if the required library is found and OFF otherwise.

Warning

Setting PTEMS_USE_TRIANGLE or PTEMS_USE_TETGEN to OFF, or PTEMS_TRIANGLE_INCDIR or PTEMS_TETGEN_INCDIR to a valid path automatically sets PTEMS_BUILD_TESTS and PTEMS_BUILD_EXAMPLES to OFF (see Included Libraries).

Warning

CMake caches the options provided. To reconfigure PTEMS with different options, it is necessary to clear option out of the cache using the extra -U NAME before the corresponding -D option.

Included Libraries

The following external libraries are included in the PTEMS repository at fixed versions under the terms of their licences:

Library	Version	License
Triangle	1.6	see License
TetGen	1.4.3 [TetGen]	MIT with Commercial Restriction

Note that the included versions of these libraries contain only a subset of the files provided by a full distribution (only the files necessary to build as a library), plus an additional CMakeLists.txt file for use in building with CMake.

Note

TetGen is not currently used (3D meshes not yet implemented).

Instead of using the built-in version of Triangle or TetGen externally built versions of Triangle and TetGen can be used by specifying the folder containing the include files (triangle.h or tetgen.h) using the option PTEMS_TRIANGLE_INCDIR and PTEMS_TETGEN_INCDIR, respectively. When linking the program using PTEMS ensure that the compiled objects or library for Triangle/TetGen are linked manually.

Triangle and TetGen can be disabled completely by specifying PTEMS_USE_TRIANGLE=NO or PTEMS_USE_TETGEN=NO, respectively.

Note

If Triangle or TetGen are disabled certain mesh construction functions (such as CreateUnstructuredMesh) will be unavailable.

Warning

Setting PTEMS_USE_TRIANGLE or PTEMS_USE_TETGEN to OFF, or PTEMS_TRIANGLE_INCDIR or PTEMS_TETGEN_INCDIR to a valid path automatically sets PTEMS_BUILD_TESTS and PTEMS_BUILD_EXAMPLES to OFF as they cannot automatically find the binary objects.

Usage

The usage of PTEMS is currently fairly low level as a C++ library. High-level interfaces are in progress. The following will demonstrate how to use PTEMS for a simple Poisson problem both a discontinuous Galerkin (DG) finite element method and a conforming virtual element method (VEM). More detailed examples can be found in the examples directory.

As a C++ library, the program using PTEMS should be set up to find the PTEMS header and built shared library files.

Headers

In the source file, first include the key headers files

#include "ptems/discretefunction.hpp"
#include "ptems/mesh.hpp"
#include "ptems/quadrature.hpp"

and a header file for the discrete space type; e.g,

#include "ptems/spaces/discontinuousfem.hpp"
#include "ptems/spaces/virtualelement.hpp"

for DG or VEM, respectively.

Construct Space

Within the main function, the following should be performed.

Create domain using ptems::DomainDefinition2D; for example, to define the unit square domain \(\Omega=(0,1)^2\) use

auto domain = ptems::DomainDefinition2D::Square();

Create a parallel processing communication object (this is necessary even in single thread mode; note that MPI is still experimental):

ptems::ParallelComm comm;

Construct a mesh from the domain:

• Uniform mesh of 4x4 squares:

  auto mesh = ptems::CreateUniformMesh(4, 4, UniformElementType2D::Square); 

• Unstructured mesh of approximately 20 triangles:

  ptems::DelaunayOptions options;
  options.MaximumTriArea = 0.05;
  options.MinimumAngle = 20;
  auto mesh = ptems::CreateUnstructuredMesh(comm, domain, options);

• Voronoi mesh of approximately 16 Voronoi:

  mesh = ptems::CreateVoronoiMesh(comm, domain, 16);

Set initial order of approximation/polynomial degree on each element:

auto polydeg = 2;

Create functions/lambdas to define the right-hand side forcing functions and Dirichlet boundary conditions (in this case, simply zero):

auto f = [](const ptems::Vector<2>& pt) { return 0; }
auto ud = [](const ptems::Vector<2>& pt) { return 0; }

Construct a DG space

auto Vh = ptems::MakeDGSpace(mesh, polydeg);

or a VEM space

auto Vh = ptems::MakeH1ConformingVESpace(mesh, 2, ud)

Note, the VEM space requires the boundary condition at construction, as opposed to DG, as it strongly imposes boundary conditions

Main loop

We then iterate for a prescribed number of iterations solving the method on a sequence of uniformly refined meshes using the following:

for(std::size_t i = 0; i != iterations; ++i)
{
    auto uh = SolveSystem(comm, mesh, Vh, f, ud);
 
    mesh->UniformRefine();
}

Here, SolveSystem is a function we need to define, see Solve System.

Within the loop body we can also compute errors, a posteriori error bounds, (see examples for details), and output the solution to VTK files:

std::ostringstream filename;
filename << "mesh"<< i << ".vtk";
mesh->Save(comm, filename.str(), ptems::MeshFileFormat::VTK_Legacy, "polydeg",
    Vh->PolynomialDegreeFunction(), "uh", uh);

Alternatively, instead of uniform mesh refinement we can perform uniform polynomial refinement:

++polydeg;
Vh->UniformPolynomialDegree(polydeg);

If we have a variable errIndicators of type std::vector<double> (or any STL compatible sequence) containing the a posteriori error indicators we can perform h- or p-adaptive mesh refinement:

mesh->Adapt(ptems::FixedFractionMarking{errIndicators, 0.25, 0.1});
Vh->Adapt(ptems::FixedFractionMarking{errIndicators, 0.25, 0.1}, polydeg);

Or even hp-adaptive:

auto hpRefinement = ptems::CreateHPPredictedErrorStrategy(Vh, 4, 0.4, 1, polydeg);
hpRefinement->AdaptWithStrategyAndDefaultSmoothing<ptems::FixedFractionMarking>(errIndicators, uh, 0.25);

Note, that the hpRefinement variable should be created once outside the loop, and additionally the following headers may be needed for adaptive refinement:

#include "ptems/refinements.hpp"
#include "ptems/hpstrategies.hpp"

Solve System

We have deferred until now the discussion of the main important routine SolveSystem. This contains the bulk of the program, and should eventually be replaced with high-level alternatives (although, the low level option will remain).

The general signature of the function is as follows:

template<typename Space, typename F, typename G>
typename ptems::DiscreteFunctionType<Space>::type SolveSystem(ptems::ParallelComm& comm,
    const ptems::PFEMesh<2>& mesh, const Space& Vh, F f, G ud)

Note, for simplicity, we pass the Space type (DG or VEM) and the boundary condition as templated arguments so we don't have to know the exact type. From the space type we can deduce the type of the discrete function (which is returned) using the ptems::DiscreteFunctionType helper struct:

typename ptems::DiscreteFunctionType<Space>::type

Similarly, we can deduce it from the variable Vh by:

typename ptems::DiscreteFunctionType<decltype(Vh)>::type

Within the body of the code, we first set initialise an empty variable of the specified type (set to zero):

typename ptems::DiscreteFunctionType<Space>::type uh(Vh, 0);

Warning

It is important to pass the 0 argument for VEM (even though it is an optional argument), as this forces the constructed function to impose any strongly imposed boundary conditions (not passing 0 does not impose the boundary conditions, which may be desirable for iterative nonlinear solvers).

It is then necessary to create the sparsity pattern for the matrix used by the solver. This is the first candidate for a function to move into PTEMS. Note, that the sparsity pattern is different depending on the method type.

• Conforming VEM with strongly imposed boundary conditions: It is first necessary to compute the real global degree of freedom index for each degree of freedom (and whether the degree of freedom is removed for strongly imposed boundary conditions):

  std::size_t sizeFullSystem = Vh->NumberGlobalDoFs();
  std::vector<std::size_t> globalDofMap(sizeFullSystem);
  auto sizeSolver = 0;
  auto strongBCs = Vh->StronglyImposedDoFs();
  auto strongBCsEnd = strongBCs.end();
  for(std::size_t i = 0; i != globalDofMap.size(); ++i)
  {
      auto it = strongBCs.find(i);
      if(it != strongBCsEnd)
      {
          globalDofMap[i] = sizeFullSystem;
      }
      else
      {
          globalDofMap[i] = sizeSolver;
          ++sizeSolver;
      }
  }

  After this the variable globalDofMap will map the index of each global degree of freedom into an index in the actual matrix, or to sizeFullSystem to denote the degree of freedom was removed as it represents a strongly imposed boundary condition degree of freedom. It is now possible to construct the sparsity pattern:

  ptems::SparsityPattern sparsity;
  for(std::size_t ele = 0; ele != mesh->ElementCount(); ++ele)
  {
      auto dofMapping = Vh->DoFMapping(ele);
      auto end = dofMapping.end();
      for(auto dof = dofMapping.begin(); dof != end; ++dof)
      {
          auto actualDoF = globalDofMap[*dof];
          if(actualDoF < sizeFullSystem)
          {
              sparsity.emplace(actualDoF, actualDoF);
              for(auto dof2 = std::next(dof); dof2 != end; ++dof2)
              {
                  auto actualDoF2 = globalDofMap[*dof2];
                  if(actualDoF2 < sizeFullSystem)
                  {
                      sparsity.emplace(actualDoF, actualDoF2);
                      sparsity.emplace(actualDoF2, actualDoF);
                  }
              }
          }
      }
  }

• Discontinuous Galerkin Finite Element Methods: Boundary conditions are now weakly imposed, so all degrees of freedom are kept; however, we now need to worry about interactions between degrees of freedom from neighbouring elements:

  std::size_t sizeSolver = Vh->NumberGlobalDoFs();
  ptems::SparsityPattern sparsity;
  for(std::size_t ele = 0; ele != mesh->ElementCount(); ++ele)
  {
      auto dofMapping = Vh->DoFMapping(ele);
      auto end = dofMapping.end();
      for(auto dof = dofMapping.begin(); dof != end; ++dof)
      {
          sparsity.emplace(*dof, *dof);
          for(auto dof2 = std::next(dof); dof2 != end; ++dof2)
          {
              sparsity.emplace(*dof, *dof2);
              sparsity.emplace(*dof2, *dof);
          }
      }
      auto element = mesh->Element(ele);
      for(std::size_t i = 0; i != element->FaceCount(); ++i)
      {
          auto right = element->Face(i)->Right();
          if(right)
          {
              auto dofMapping2 = Vh->DoFMapping(right->Index());
              auto end2 = dofMapping2.end();
              for(auto dof = dofMapping.begin(); dof != end; ++dof)
              {
                  for(auto dof2 = dofMapping2.begin(); dof2 != end2; ++dof2)
                  {
                      sparsity.emplace(*dof, *dof2);
                      sparsity.emplace(*dof2, *dof);
                  }
              }
          }
      }
  }

Now we have the sparsity pattern, we can construct the linear solver and right-hand side vector:

ptems::MumpsSolver<double> solver(comm, sparsity, sizeSolver, ptems::MumpsMatrixType::Unsymmetric, false);
std::vector<double> rhs(sizeSolver);

The last argument is designed for use for distributed construction when using MPI (should be ignored for single process, but we set to false anyway).

We now have a matrix and right-hand side containing zeros, so we need to iterate over the elements and populate the matrix and right-hand side using numerical integration. In order to do this we first initialise a cache object for caching the quadrature points on the elements and weights, as PTEMS dynamically generates these:

ptems::GaussLegendreQuadraturePoints<2> quadElements;

Additionally, for DG we need edge integration, so we need a cache for edge quadrature points:

ptems::GaussLegendreQuadraturePoints<1> quadFaces;

We now loop through each element and perform numerical integration of \((\nabla \phi_j, \nabla \phi_i)_K\) and \((f,\phi_i)_K\) for each pair of degrees of freedom of index i and j. The following is the code for VEM:

for(std::size_t ele = 0; ele != mesh->ElementCount(); ++ele)
{
    auto element = mesh->Element(ele);
    auto p = Vh->PolynomialDegree(ele);
    auto quadPts = element->ComputeQuadraturePoints(quadElements, p[0]+4);
    auto dofMapping = Vh->DoFMapping(ele);
    for(std::size_t k = 0; k != quadPts.size(); ++k)
    {
        auto globalPt = element->LocalToGlobal(quadPts[k].Point);
        auto basis = Vh->GradBasisFunctions(ele, quadPts[k].Point);
        for(std::size_t i = 0; i != basis.size(0); ++i)
        {
            auto dofI = globalDofMap[dofMapping(0,i)];
            if(dofI < sizeFullSystem)
            {
                for(std::size_t j = 0; j != basis.size(0); ++j)
                {
                    auto dofJ = globalDofMap[dofMapping(0,j)];
                    if(dofJ < sizeFullSystem)
                    {
                        solver(dofI,dofJ) += quadPts[k].Weight*
                            basis(0,i).Gradient.Dot(basis(0,j).Gradient);
                    }
                    else
                    {
                        // Boundary DoF - strongly impose in RHS
                        rhs[dofI] -= quadPts[k].Weight*uh[dofMapping(0,j)]
                            *basis(0,i).Gradient.Dot(basis(0,j).Gradient);
                    }
                }
                rhs[dofI] += quadPts[k].Weight*f(globalPt)*basis(0,i).Func;
            }
        }
    }
}

For DG, we simply replace the lines

auto dofI = globalDofMap[dofMapping(0,i)];
if(dofI < sizeFullSystem)
{
    ...
}

with

auto dofI = dofMapping(0,i);
...

and similarly for dofJ; additionally, remove the strongly imposed boundary conditions case

rhs[dofI] -= quadPts[k].Weight*uh[dofMapping(0,j)]
    *basis(0,i).Gradient.Dot(basis(0,j).Gradient);

This constructs the matrix and right-hand side for the volume integrations, but both methods require additional terms:

• Virtual element method uses a stabilisation, in this case we consider the traditional dofi-dofi stabilisation. Inside the loop over elements we add the following:

  auto dofOfBasis = Vh->EvaluateDoFAtBasis(ele);
  auto dofOfBasisProj = Vh->EvaluateDoFAtBasisProjection(ele);
  for(std::size_t k = 0; k != dofOfBasis.size(0); ++k)
  {
      const auto& dofBasisK = dofOfBasis(0,k);
      const auto& dofBasisProjK = dofOfBasisProj(0,k);
      for(std::size_t i = 0; i != dofBasisK.size(); ++i)
      {
          auto dofI = globalDofMap[dofMapping(0,i)];
          if(dofI < sizeFullSystem)
          {
              for(std::size_t j = 0; j != dofBasisK.size(); ++j)
              {
                  auto dofJ = globalDofMap[dofMapping(0,j)];
                  if(dofJ < sizeFullSystem)
                  {
                      solver(dofI,dofJ) += (dofBasisK[i]-dofBasisProjK[i])
                          * (dofBasisK[j]-dofBasisProjK[j]);
                  }
              }
          }
      }
  }

• Discontinuous Galerkin finite element method (in this case interior penalty, where we have some parameter theta set to -1, 0, or 1 for symmetric, incomplete, and nonsymmetric, respectively), requires face integral terms. Inside the loop over elements we add the following, where interiorPenalty is set to the constant to use for the interior penalty term:

  for(std::size_t e = 0; e != element->FaceCount(); ++e)
  {
      auto face = element->Face(e);
      auto quadPtsEdge = face->ComputeQuadraturePoints(quadFaces, p[0]+4);
      auto right = face->Right();
      auto n = face->Normal();
      std::optional<decltype(Vh->DoFMapping(ele))> dofMapping2;
      double pf = static_cast<double>(p[0]);
      double edgeAvg = 1.0;
      double hf = element->Volume()/face->Volume();
      if(right)
      {
          dofMapping2 = Vh->DoFMapping(right->Index());
          pf = static_cast<double>(std::max(p[0],Vh->PolynomialDegree(right->Index())[0]));
          edgeAvg = 0.5;
          hf = std::min(hf, right->Volume()/face->Volume());
      }
      for(std::size_t k = 0; k != quadPtsEdge.size(); ++k)
      {
          auto globalPt = face->LocalToGlobal(quadPtsEdge[k].Point);
          auto localPt1 = element->GlobalToLocal(globalPt);
          auto basis1 = Vh->GradBasisFunctions(ele, localPt1);
          for(std::size_t i = 0; i != basis1.size(0); ++i)
          {
              // Interior: +/+ terms
              for(std::size_t j = 0; j != basis1.size(0); ++j)
              {
                  solver(dofMapping(0,i),dofMapping(0,j)) += quadPtsEdge[k].Weight*(
                      -edgeAvg*(basis1(0,j).Gradient*basis1(0,i).Func
                          -theta*basis1(0,i).Gradient*basis1(0,j).Func).Dot(n)
                      +penaltyParam*basis1(0,i).Func*basis1(0,j).Func*pf*pf/hf
                  );
              }
          }
          if(right)
          {
              // Interior: -/+ terms
              auto localPt2 = right->GlobalToLocal(globalPt);
              auto basis2 = Vh->GradBasisFunctions(right->Index(), localPt2);
              for(std::size_t i = 0; i != basis1.size(0); ++i)
              {
                  for(std::size_t j = 0; j != basis2.size(0); ++j)
                  {
                      solver(dofMapping(0,i),(*dofMapping2)(0,j)) += quadPtsEdge[k].Weight*(
                          -edgeAvg*(basis2(0,j).Gradient*basis1(0,i).Func
                              +theta*basis1(0,i).Gradient*basis2(0,j).Func).Dot(n)
                          -penaltyParam*basis1(0,i).Func*basis2(0,j).Func*pf*pf/hf
                      );
                  }
              }
          }
          else
          {
              auto udx = ud(globalPt);
              // Boundary
              for(std::size_t i = 0; i != basis1.size(0); ++i)
              {
                  rhs[dofMapping(0,i)] += quadPtsEdge[k].Weight*udx*(
                      theta*basis1(0,i).Gradient.Dot(n)
                      +penaltyParam*basis1(0,i).Func*pf*pf/hf
                  );
              }
          }
      }
  }

  Note, for average and jumps, as we are doing this for every element we consider every interior edge twice. Therefore, we only consider the so-called plus-plus and plus-minus terms; i.e., the trial function from the inside of the element but the test from both the element and its neighbour.

Now, we just need to solve and return the result. This is slightly different for the two methods again, due to the fact that DG gives the full solution (no strongly imposed boundaries); whereas, VEM doesn't give the boundary dofs:

• Strongly imposed boundaries: We need to compute the result vector and then copy into the solution function at the correct indices as the solution already contains the correct boundary values:

  auto result = solver.Solve(rhs);
  for(std::size_t i = 0; i != globalDofMap.size(); ++i)
  {
      if(globalDofMap[i] < sizeFullSystem)
      {
          uh[i] = result[globalDofMap[i]];
      }
  }
  return uh;

• Weakly imposed boundaries: Here we can directly assign the result to the solution function:

  return (uh = solver.Solve(rhs));

Nonlinear Solver

The above is for a linear problem. For a nonlinear problem a damped Newton solver is provided, which requires additional steps.

Note

The nonlinear solver is designed to implement a damped Newton solver; however, it can easily be used for any damped or undamped Newton-like solver of the following form:

\[u_{n+1} = u_{n} - \varepsilon A^{-1}R(u_n)\]

where \(\varepsilon\in(0,1]\) is the damping parameter, which the algorithm may alter, R is the residual operator used to detect convergence and A is a linear operator.

To solve a nonlinear system we change the SolveSystem implementation. The key ingredients are as follows:

• Create a nonlinear solver, in addition to the linear solver:

  ptems::DampingParameters<double> dampingParams{0.2, 0.2, 0.01};
  ptems::NonlinearSolver<double> nlSolver(dampingParams);

  Note, the three arguments to the DampingParameters constructor are the initial damping, increment, and minimum; therefore passing zero increment and initial == minimum removes damping; alternatively, for no damping and damping parameter of 1, simply use the default constructor to NonlinearSolver.
• Create a function (lambda function is easiest) to compute the residual:

  auto constructResidual = [&](typename ptems::DiscreteFunctionType<Space>::type& phi,
          std::vector<SystemType>& residual)

  This function takes the solution from the previous step (phi), and a vector (of length of the system size minus strongly imposed boundary) to populate with the residual vector. The contents of this function is similar to the linear solver case; except, only one basis (test) functions are needed.

• Create a function to solve the linear system \(\langle A(u),v\rangle = \langle R(\phi), v \rangle\) (e.g., Jacobian for Newton):

  auto solveJacobian = [&](const typename ptems::DiscreteFunctionType<Space>::type& phi,
          const std::vector<SystemType>& residual)

  This function takes the solution from the previous step (phi), and the residual returned by the previous function, and should return the solution of the linearised system at this step. This function is, therefore, almost identical to the linear case.

  Warning

  When constructing the solution for the linear step it is important not to pass the second argument to the constructor in the case of strongly imposed boundary conditions as the update must have zero boundary conditions:

  typename ptems::DiscreteFunctionType<Space>::type uh(Vh);

• Then, we pass these two functions to the nonlinear solver, along with the initial guess which we enforce the strongly imposed boundary conditions on in the case of VEM:

  typename ptems::DiscreteFunctionType<Space>::type uh(Vh, 0);
  return nlSolver.Solve(solveJacobian, constructResidual, uh, 1e-8, 100);

  The last two optional arguments specify the desired tolerance in the residual and the maximum number of iterations. The function SolveWithStatus can be used instead where the tolerance and the number of iterations are passed by mutable reference, and are updated with the achieved tolerance and actual number of iterations, respectively.

Current Issues

Below are listed some known issues and experimental features:

• The ability to combine multiple discrete spaces into a Cartesian Product Space (ptems::DiscreteCartesianProductSpace) is broken and should not be used (issue #1)
• 64-bit builds on Windows are broken due to issues with Triangle (issue #2) - an experimental fix is available in the win64 branch
• Examples with voronoi meshes occasionally fail (issue #3)
• MPI support for multi-processors is extremely patchy and should not be used
• Adaptivity for agglomerated elements is broken and not support yet
• Some inconsistent crashes and bugs seem to occur when the mesh is very heavily refined, likely caused by rounding errors.
• GMRES linear solver is extremely experimental

Roadmap

The following describes the list of future features, in decreasing order of importance, to be implemented:

• Support for space-time spaces (experimental spacetime branch exists)
• Add r-tree method for mesh agglomeration and refinement (https://arxiv.org/abs/2404.18505)
• Support for 3D meshes; including 3D Voronoi construction and element agglomeration
• Add ability to import/export different mesh formats
• Add support for mesh refinement with agglomerated elements
• Symbolic (or similar) method to construct the system to solve
• MPI support
• Python wrapper

License Information

The PTEMS project is licensed under the BSD 3-Clause license, see License.

The libraries used by PTEMS may use other licenses. Notably, two libraries are included with the PTEMS source (Triangle and TetGen), which contain licenses which are ONLY valid for private, research, or educational use; i.e., NOT FOR COMMERCIAL USE (unless permission is obtained from the respective copyright owners). See LICENSE for more details.

Note, that Triangle and TetGen can be removed from PTEMS (or replaced with different versions), see Included Libraries.

Authors

• Scott Congreve (congr.nosp@m.eve@.nosp@m.karli.nosp@m.n.mf.nosp@m.f.cun.nosp@m.i.cz)

---

AppleClang

Apple Clang 11.0.0 (available in XCode 11 or later) on macOS.

MSVC

Visual Studio support is currently untested

List

A CMake list is a list of semicolon ';' separated strings. Alternatively the option can be specified multiple times for each item in the list.

TetGen

Note this is an old version of Tetgen (from 2011). This is included as it is available under a (no commerical use without permission) version of the MIT license which is mostly compatible with PTEMS' BSD-3 license (except the no commercial restriction); whereas, up-to-date versions appear to use GNU Affero GPL, or a paid commerical license, which is NOT compatible.