No Matches
Classes | Functions | Variables
analytical_poisson.cpp File Reference
#include <BasicFiniteElements.hpp>
#include <PoissonOperators.hpp>
#include <AuxPoissonFunctions.hpp>

Go to the source code of this file.


struct  ExactFunction
 Function. More...
struct  ExactFunctionGrad
 Exact gradient. More...
struct  ExactLaplacianFunction
 Laplacian of function. More...


int main (int argc, char *argv[])


static char help [] = "...\n\n"

Function Documentation

◆ main()

int main ( int  argc,
char *  argv[] 

< Volume element for the matrix

< Boundary element for the matrix

< Volume element to assemble vector

< Volume element to assemble vector

< Volume element evaluate error

< Volume element to Post-process results

< Null element do nothing

Definition at line 70 of file analytical_poisson.cpp.

70 {
72 // Initialize PETSc
73 MoFEM::Core::Initialize(&argc,&argv,(char *)0,help);
75 try {
77 // Create MoAB database
78 moab::Core moab_core; // create database
79 moab::Interface &moab = moab_core; // create interface to database
81 // Get command line options
82 int order = 3; // default approximation order
83 PetscBool flg_test = PETSC_FALSE; // true check if error is numerical error
84 CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "Poisson's problem options",
85 "none");
86 // Set approximation order
87 CHKERR PetscOptionsInt("-order", "approximation order", "", order, &order,
89 // Set testing (used by CTest)
90 CHKERR PetscOptionsBool("-test", "if true is ctest", "", flg_test,
91 &flg_test, PETSC_NULL);
92 ierr = PetscOptionsEnd();
95 // Create MoFEM database and link it to MoAB
96 MoFEM::Core mofem_core(moab); // create database
97 MoFEM::Interface &m_field = mofem_core; // create interface to database
98 // Register DM Manager
99 CHKERR DMRegister_MoFEM("DMMOFEM"); // register MoFEM DM in PETSc
101 // Create vector to store approximation global error
102 Vec global_error;
103 CHKERR PoissonExample::AuxFunctions(m_field).createGhostVec(&global_error);
105 // First we crate elements, implementation of elements is problem independent,
106 // we do not know yet what fields are present in the problem, or
107 // even we do not decided yet what approximation base or spaces we
108 // are going to use. Implementation of element is free from
109 // those constrains and can be used in different context.
111 // Elements used by KSP & DM to assemble system of equations
112 boost::shared_ptr<ForcesAndSourcesCore> domain_lhs_fe; ///< Volume element for the matrix
113 boost::shared_ptr<ForcesAndSourcesCore> boundary_lhs_fe; ///< Boundary element for the matrix
114 boost::shared_ptr<ForcesAndSourcesCore> domain_rhs_fe; ///< Volume element to assemble vector
115 boost::shared_ptr<ForcesAndSourcesCore> boundary_rhs_fe; ///< Volume element to assemble vector
116 boost::shared_ptr<ForcesAndSourcesCore> domain_error; ///< Volume element evaluate error
117 boost::shared_ptr<PoissonExample::PostProcFE>
118 post_proc_volume; ///< Volume element to Post-process results
119 boost::shared_ptr<ForcesAndSourcesCore> null; ///< Null element do nothing
120 {
121 // Add problem specific operators the generic finite elements to calculate matrices and vectors.
124 ExactFunction(), ExactLaplacianFunction(), domain_lhs_fe,
125 boundary_lhs_fe, domain_rhs_fe, boundary_rhs_fe);
126 // Add problem specific operators the generic finite elements to calculate error on elements and global error
127 // in H1 norm
130 global_error, domain_error);
131 // Post-process results
133 .creatFEToPostProcessResults(post_proc_volume);
134 }
136 // Get simple interface is simplified version enabling quick and
137 // easy construction of problem.
138 Simple *simple_interface;
139 // Query interface and get pointer to Simple interface
140 CHKERR m_field.getInterface(simple_interface);
142 // Build problem with simple interface
143 {
145 // Get options for simple interface from command line
146 CHKERR simple_interface->getOptions();
147 // Load mesh file to database
148 CHKERR simple_interface->loadFile();
150 // Add field on domain and boundary. Field is declared by space and base and rank. space
151 // can be H1. Hcurl, Hdiv and L2, base can be AINSWORTH_LEGENDRE_BASE, DEMKOWICZ_JACOBI_BASE and more,
152 // where rank is number of coefficients for dof.
153 //
154 // Simple interface assumes that there are four types of field; 1) defined
155 // on body domain, 2) fields defined on body boundary, 3) skeleton field defined
156 // on finite element skeleton. Finally data field is defined on body domain. Data field
157 // is not solved but used for post-process or to keep material parameters, etc. Here
158 // we using it to calculate approximation error on elements.
160 // Add domain filed "U" in space H1 and Legendre base, Ainsworth recipe is used
161 // to construct base functions.
162 CHKERR simple_interface->addDomainField("U",H1,AINSWORTH_LEGENDRE_BASE,1);
163 // Add Lagrange multiplier field on body boundary
164 CHKERR simple_interface->addBoundaryField("L",H1,AINSWORTH_LEGENDRE_BASE,1);
165 // Add error (data) field, we need only L2 norm. Later order is set to 0, so this
166 // is piecewise discontinuous constant approx., i.e. 1 DOF for element. You can use
167 // more DOFs and collate moments of error to drive anisotropic h/p-adaptivity, however
168 // this is beyond this example.
169 CHKERR simple_interface->addDataField("ERROR",L2,AINSWORTH_LEGENDRE_BASE,1);
171 // Set fields order domain and boundary fields.
172 CHKERR simple_interface->setFieldOrder("U",order); // to approximate function
173 CHKERR simple_interface->setFieldOrder("L",order); // to Lagrange multipliers
174 CHKERR simple_interface->setFieldOrder("ERROR",0); // approximation order for error
176 // Setup problem. At that point database is constructed, i.e. fields, finite elements and
177 // problem data structures with local and global indexing.
178 CHKERR simple_interface->setUp();
180 }
182 // Get access to PETSC-MoFEM DM manager.
183 // or more derails see <http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/DM/index.html>
184 // Form that point internal MoFEM data structures are linked with PETSc interface. If
185 // DM functions contains string MoFEM is is MoFEM specific DM interface function,
186 // otherwise DM function part of standard PETSc interface.
188 DM dm;
189 // Get dm
190 CHKERR simple_interface->getDM(&dm);
192 // Set KSP context for DM. At that point only elements are added to DM operators.
193 // Calculations of matrices and vectors is executed by KSP solver. This part
194 // of the code makes connection between implementation of finite elements and
195 // data operators with finite element declarations in DM manager. From that
196 // point DM takes responsibility for executing elements, calculations of
197 // matrices and vectors, and solution of the problem.
198 {
199 // Set operators for KSP solver
201 dm, simple_interface->getDomainFEName(), domain_lhs_fe, null, null);
203 dm, simple_interface->getBoundaryFEName(), boundary_lhs_fe, null,
204 null);
205 // Set calculation of the right hand side vector for KSP solver
206 CHKERR DMMoFEMKSPSetComputeRHS(dm, simple_interface->getDomainFEName(),
207 domain_rhs_fe, null, null);
208 CHKERR DMMoFEMKSPSetComputeRHS(dm, simple_interface->getBoundaryFEName(),
209 boundary_rhs_fe, null, null);
210 }
212 // Solve problem, only PETEc interface here.
213 {
215 // Create the right hand side vector and vector of unknowns
216 Vec F,D;
217 CHKERR DMCreateGlobalVector(dm,&F);
218 // Create unknown vector by creating duplicate copy of F vector. only
219 // structure is duplicated no values.
220 CHKERR VecDuplicate(F,&D);
222 // Create solver and link it to DM
223 KSP solver;
225 CHKERR KSPSetFromOptions(solver);
226 CHKERR KSPSetDM(solver,dm);
227 // Set-up solver, is type of solver and pre-conditioners
228 CHKERR KSPSetUp(solver);
229 // At solution process, KSP solver using DM creates matrices, Calculate
230 // values of the left hand side and the right hand side vector. then
231 // solves system of equations. Results are stored in vector D.
232 CHKERR KSPSolve(solver,F,D);
234 // Scatter solution on the mesh. Stores unknown vector on field on the mesh.
237 // Clean data. Solver and vector are not needed any more.
238 CHKERR KSPDestroy(&solver);
239 CHKERR VecDestroy(&D);
240 CHKERR VecDestroy(&F);
241 }
243 // Calculate error
244 {
245 // Loop over all elements in mesh, and run users operators on each element.
246 CHKERR DMoFEMLoopFiniteElements(dm, simple_interface->getDomainFEName(),
247 domain_error);
249 global_error);
250 CHKERR PoissonExample::AuxFunctions(m_field).printError(global_error);
251 if (flg_test == PETSC_TRUE) {
252 CHKERR PoissonExample::AuxFunctions(m_field).testError(global_error);
253 }
254 }
256 {
257 // Loop over all elements in the mesh and for each execute post_proc_volume
258 // element and operators on it.
259 CHKERR DMoFEMLoopFiniteElements(dm, simple_interface->getDomainFEName(),
260 post_proc_volume);
261 // Write results
262 post_proc_volume->writeFile("out_vol.h5m");
263 }
265 // Destroy DM, no longer needed.
266 CHKERR DMDestroy(&dm);
268 // Destroy ghost vector
269 CHKERR VecDestroy(&global_error);
270 }
273 // finish work cleaning memory, getting statistics, etc.
276 return 0;
static char help[]
Catch errors.
Definition: definitions.h:372
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
@ L2
field with C-1 continuity
Definition: definitions.h:88
@ H1
continuous field
Definition: definitions.h:85
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:483
#define CHKERR
Inline error check.
Definition: definitions.h:535
constexpr int order
@ F
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:47
PetscErrorCode DMMoFEMKSPSetComputeRHS(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set KSP right hand side evaluation function
Definition: DMMoFEM.cpp:623
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition: DMMoFEM.cpp:572
PetscErrorCode DMoFEMMeshToGlobalVector(DM dm, Vec g, InsertMode mode, ScatterMode scatter_mode)
set ghosted vector values on all existing mesh entities
Definition: DMMoFEM.cpp:521
PetscErrorCode DMMoFEMKSPSetComputeOperators(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
Set KSP operators and push mofem finite element methods.
Definition: DMMoFEM.cpp:664
double D
const FTensor::Tensor2< T, Dim, Dim > Vec
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
Exact gradient.
Core (interface) class.
Definition: Core.hpp:82
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:72
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
Deprecated interface functions.
Simple interface for fast problem set-up.
Definition: Simple.hpp:27
const std::string getBoundaryFEName() const
Get the Boundary FE Name.
Definition: Simple.hpp:348
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name, LoadFileFunc loadFunc=defaultLoadFileFunc)
Load mesh file.
Definition: Simple.cpp:194
MoFEMErrorCode addDataField(const std::string &name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_ZERO, int verb=-1)
Add field on domain.
Definition: Simple.cpp:320
MoFEMErrorCode addDomainField(const std::string &name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_ZERO, int verb=-1)
Add field on domain.
Definition: Simple.cpp:264
MoFEMErrorCode getOptions()
get options
Definition: Simple.cpp:180
MoFEMErrorCode getDM(DM *dm)
Get DM.
Definition: Simple.cpp:667
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
Definition: Simple.cpp:473
MoFEMErrorCode addBoundaryField(const std::string &name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_ZERO, int verb=-1)
Add field on boundary.
Definition: Simple.cpp:282
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
Definition: Simple.cpp:611
const std::string getDomainFEName() const
Get the Domain FE Name.
Definition: Simple.hpp:341
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
MoFEMErrorCode createGhostVec(Vec *ghost_vec) const
MoFEMErrorCode testError(Vec ghost_vec)
Test error.
MoFEMErrorCode assembleGhostVector(Vec ghost_vec) const
Assemble error vector.
MoFEMErrorCode printError(Vec ghost_vec)
Print error.
Create finite elements instances.
MoFEMErrorCode creatFEToPostProcessResults(boost::shared_ptr< PostProcFE > &post_proc_volume) const
Create finite element to post-process results.
MoFEMErrorCode createFEToEvaluateError(boost::function< double(const double, const double, const double)> f_u, boost::function< FTensor::Tensor1< double, 3 >(const double, const double, const double)> g_u, Vec global_error, boost::shared_ptr< ForcesAndSourcesCore > &domain_error) const
Create finite element to calculate error.
MoFEMErrorCode createFEToAssembleMatrixAndVector(boost::function< double(const double, const double, const double)> f_u, boost::function< double(const double, const double, const double)> f_source, boost::shared_ptr< ForcesAndSourcesCore > &domain_lhs_fe, boost::shared_ptr< ForcesAndSourcesCore > &boundary_lhs_fe, boost::shared_ptr< ForcesAndSourcesCore > &domain_rhs_fe, boost::shared_ptr< ForcesAndSourcesCore > &boundary_rhs_fe, bool trans=true) const
Create finite element to calculate matrix and vectors.

Variable Documentation

◆ help

char help[] = "...\n\n"

Definition at line 18 of file analytical_poisson.cpp.