v0.9.0
Classes | Functions | Variables
analytical_nonlinear_poisson.cpp File Reference
#include <BasicFiniteElements.hpp>
#include <PoissonOperators.hpp>
#include <AuxPoissonFunctions.hpp>

Go to the source code of this file.

Classes

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

Functions

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

Variables

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

Function Documentation

◆ main()

int main ( int  argc,
char *  argv[] 
)
Examples
analytical_nonlinear_poisson.cpp.

Definition at line 98 of file analytical_nonlinear_poisson.cpp.

98  {
99 
100  // Initialize PETSc
101  MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
102  // Create MoAB database
103  moab::Core moab_core; // create database
104  moab::Interface &moab = moab_core; // create interface to database
105 
106  try {
107 
108  // Get command line options
109  int order = 3; // default approximation order
110  PetscBool flg_test = PETSC_FALSE; // true check if error is numerical error
111  CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "Poisson's problem options",
112  "none");
113  // Set approximation order
114  CHKERR PetscOptionsInt("-order", "approximation order", "", order, &order,
115  PETSC_NULL);
116  // Set testing (used by CTest)
117  CHKERR PetscOptionsBool("-test", "if true is ctest", "", flg_test,
118  &flg_test, PETSC_NULL);
119  ierr = PetscOptionsEnd();
120  CHKERRG(ierr)
121 
122  // Create MoFEM database and link it to MoAB
123  MoFEM::Core mofem_core(moab); // create database
124  MoFEM::Interface &m_field = mofem_core; // create interface to database
125  // Register DM Manager
126  CHKERR DMRegister_MoFEM("DMMOFEM"); // register MoFEM DM in PETSc
127 
128  // Create vector to store approximation global error
129  Vec global_error;
130  CHKERR PoissonExample::AuxFunctions(m_field).createGhostVec(&global_error);
131 
132  // First we crate elements, implementation of elements is problem
133  // independent, we do not know yet what fields are present in the problem,
134  // or even we do not decided yet what approximation base or spaces we are
135  // going to use. Implementation of element is free from those constrains and
136  // can be used in different context.
137 
138  // Elements used by KSP & DM to assemble system of equations
139  boost::shared_ptr<ForcesAndSourcesCore>
140  domain_lhs_fe; ///< Volume element for the matrix
141  boost::shared_ptr<ForcesAndSourcesCore>
142  boundary_lhs_fe; ///< Boundary element for the matrix
143  boost::shared_ptr<ForcesAndSourcesCore>
144  domain_rhs_fe; ///< Volume element to assemble vector
145  boost::shared_ptr<ForcesAndSourcesCore>
146  boundary_rhs_fe; ///< Volume element to assemble vector
147  boost::shared_ptr<ForcesAndSourcesCore>
148  domain_error; ///< Volume element evaluate error
149  boost::shared_ptr<ForcesAndSourcesCore>
150  post_proc_volume; ///< Volume element to Post-process results
151  boost::shared_ptr<ForcesAndSourcesCore> null; ///< Null element do nothing
152  {
153  // Add problem specific operators the generic finite elements to calculate
154  // matrices and vectors.
158  domain_lhs_fe, boundary_lhs_fe, domain_rhs_fe, boundary_rhs_fe,
159  VolRuleNonlinear());
160  // Add problem specific operators the generic finite elements to calculate
161  // error on elements and global error in H1 norm
164  global_error, domain_error);
165  // Post-process results
167  .creatFEToPostProcessResults(post_proc_volume);
168  }
169 
170  // Get simple interface is simplified version enabling quick and
171  // easy construction of problem.
172  Simple *simple_interface;
173  // Query interface and get pointer to Simple interface
174  CHKERR m_field.getInterface(simple_interface);
175 
176  // Build problem with simple interface
177  {
178 
179  // Get options for simple interface from command line
180  CHKERR simple_interface->getOptions();
181  // Load mesh file to database
182  CHKERR simple_interface->loadFile();
183 
184  // Add field on domain and boundary. Field is declared by space and base
185  // and rank. space can be H1. Hcurl, Hdiv and L2, base can be
186  // AINSWORTH_LEGENDRE_BASE, DEMKOWICZ_JACOBI_BASE and more, where rank is
187  // number of coefficients for dof.
188  //
189  // Simple interface assumes that there are four types of field; 1) defined
190  // on body domain, 2) fields defined on body boundary, 3) skeleton field
191  // defined on finite element skeleton. Finally data field is defined on
192  // body domain. Data field is not solved but used for post-process or to
193  // keep material parameters, etc. Here we using it to calculate
194  // approximation error on elements.
195 
196  // Add domain filed "U" in space H1 and Legendre base, Ainsworth recipe is
197  // used to construct base functions.
198  CHKERR simple_interface->addDomainField("U", H1, AINSWORTH_LEGENDRE_BASE,
199  1);
200  // Add Lagrange multiplier field on body boundary
201  CHKERR simple_interface->addBoundaryField("L", H1,
203  // Add error (data) field, we need only L2 norm. Later order is set to 0,
204  // so this is piecewise discontinuous constant approx., i.e. 1 DOF for
205  // element. You can use more DOFs and collate moments of error to drive
206  // anisotropic h/p-adaptivity, however this is beyond this example.
207  CHKERR simple_interface->addDataField("ERROR", L2,
209 
210  // Set fields order domain and boundary fields.
211  CHKERR simple_interface->setFieldOrder("U",
212  order); // to approximate function
213  CHKERR simple_interface->setFieldOrder("L",
214  order); // to Lagrange multipliers
215  CHKERR simple_interface->setFieldOrder(
216  "ERROR", 0); // approximation order for error
217 
218  // Setup problem. At that point database is constructed, i.e. fields,
219  // finite elements and problem data structures with local and global
220  // indexing.
221  CHKERR simple_interface->setUp();
222  }
223 
224  // Get access to PETSC-MoFEM DM manager.
225  // or more derails see
226  // <http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/DM/index.html>
227  // Form that point internal MoFEM data structures are linked with PETSc
228  // interface. If DM functions contains string MoFEM is is MoFEM specific DM
229  // interface function, otherwise DM function part of standard PETSc
230  // interface.
231 
232  DM dm;
233  // Get dm
234  CHKERR simple_interface->getDM(&dm);
235 
236  // Set KSP context for DM. At that point only elements are added to DM
237  // operators. Calculations of matrices and vectors is executed by KSP
238  // solver. This part of the code makes connection between implementation of
239  // finite elements and data operators with finite element declarations in DM
240  // manager. From that point DM takes responsibility for executing elements,
241  // calculations of matrices and vectors, and solution of the problem.
242  {
243  // Set operators for SNES solver
244  CHKERR DMMoFEMSNESSetJacobian(dm, simple_interface->getDomainFEName(),
245  domain_lhs_fe, null, null);
246  CHKERR DMMoFEMSNESSetJacobian(dm, simple_interface->getBoundaryFEName(),
247  boundary_lhs_fe, null, null);
248  // Set calculation of the right hand side vector for KSP solver
249  CHKERR DMMoFEMSNESSetFunction(dm, simple_interface->getDomainFEName(),
250  domain_rhs_fe, null, null);
251  CHKERR DMMoFEMSNESSetFunction(dm, simple_interface->getBoundaryFEName(),
252  boundary_rhs_fe, null, null);
253  }
254 
255  // Solve problem, only PETEc interface here.
256  {
257 
258  // Create the right hand side vector and vector of unknowns
259  Vec F, D;
260  CHKERR DMCreateGlobalVector(dm, &F);
261  // Create unknown vector by creating duplicate copy of F vector. only
262  // structure is duplicated no values.
263  CHKERR VecDuplicate(F, &D);
264 
265  // Create solver and link it to DM
266  SNES solver;
267  CHKERR SNESCreate(PETSC_COMM_WORLD, &solver);
268  CHKERR SNESSetFromOptions(solver);
269  CHKERR SNESSetDM(solver, dm);
270  // Set-up solver, is type of solver and pre-conditioners
271  CHKERR SNESSetUp(solver);
272  // At solution process, KSP solver using DM creates matrices, Calculate
273  // values of the left hand side and the right hand side vector. then
274  // solves system of equations. Results are stored in vector D.
275  CHKERR SNESSolve(solver, F, D);
276 
277  // Scatter solution on the mesh. Stores unknown vector on field on the
278  // mesh.
279  CHKERR DMoFEMMeshToGlobalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
280 
281  // Clean data. Solver and vector are not needed any more.
282  CHKERR SNESDestroy(&solver);
283  CHKERR VecDestroy(&D);
284  CHKERR VecDestroy(&F);
285  }
286 
287  // Calculate error
288  {
289  // Loop over all elements in mesh, and run users operators on each
290  // element.
291  CHKERR DMoFEMLoopFiniteElements(dm, simple_interface->getDomainFEName(),
292  domain_error);
294  global_error);
295  CHKERR PoissonExample::AuxFunctions(m_field).printError(global_error);
296  if (flg_test == PETSC_TRUE) {
297  CHKERR PoissonExample::AuxFunctions(m_field).testError(global_error);
298  }
299  }
300 
301  {
302  // Loop over all elements in the mesh and for each execute
303  // post_proc_volume element and operators on it.
304  CHKERR DMoFEMLoopFiniteElements(dm, simple_interface->getDomainFEName(),
305  post_proc_volume);
306  // Write results
307  CHKERR boost::static_pointer_cast<PostProcVolumeOnRefinedMesh>(
308  post_proc_volume)
309  ->writeFile("out_vol.h5m");
310  }
311 
312  // Destroy DM, no longer needed.
313  CHKERR DMDestroy(&dm);
314 
315  // Destroy ghost vector
316  CHKERR VecDestroy(&global_error);
317  }
318  CATCH_ERRORS;
319 
320  // finish work cleaning memory, getting statistics, etc.
322 
323  return 0;
324 }
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
Definition: Simple.cpp:507
PetscErrorCode DMMoFEMSNESSetFunction(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set SNES residual evaluation function
Definition: DMMMoFEM.cpp:626
const std::string getBoundaryFEName() const
Definition: Simple.hpp:222
MoFEMErrorCode getDM(DM *dm)
Get DM.
Definition: Simple.cpp:517
MoFEMErrorCode createFEToAssembleMatrixAndVectorForNonlinearProblem(boost::function< double(const double, const double, const double)> f_u, boost::function< double(const double, const double, const double)> f_source, boost::function< double(const double)> a, boost::function< double(const double)> diff_a, 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, ForcesAndSourcesCore::RuleHookFun vol_rule, ForcesAndSourcesCore::RuleHookFun face_rule=FaceRule(), bool trans=true) const
Create finite element to calculate matrix and vectors.
MoFEMErrorCode createGhostVec(Vec *ghost_vec) const
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:544
Core (interface) class.
Definition: Core.hpp:50
MoFEMErrorCode loadFile(const std::string options="PARALLEL=READ_PART;" "PARALLEL_RESOLVE_SHARED_ENTS;" "PARTITION=PARALLEL_PARTITION;")
Load mesh file.
Definition: Simple.cpp:62
static char help[]
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method)
Executes FEMethod for finite elements in DM.
Definition: DMMMoFEM.cpp:496
MoFEMErrorCode testError(Vec ghost_vec)
Test error.
Simple interface for fast problem set-up.
Definition: Simple.hpp:36
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:51
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:145
MoFEMErrorCode getInterface(const MOFEMuuid &uuid, IFACE *&iface) const
Get interface by uuid and return reference to pointer of interface.
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:116
PetscErrorCode DMoFEMMeshToGlobalVector(DM dm, Vec g, InsertMode mode, ScatterMode scatter_mode)
set ghosted vector values on all existing mesh entities
Definition: DMMMoFEM.cpp:446
Create finite elements instances.
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:53
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:146
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:101
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.
PetscErrorCode DMMoFEMSNESSetJacobian(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set SNES Jacobian evaluation function
Definition: DMMMoFEM.cpp:667
#define CHKERR
Inline error check.
Definition: definitions.h:596
const std::string getDomainFEName() const
Definition: Simple.hpp:221
MoFEMErrorCode creatFEToPostProcessResults(boost::shared_ptr< ForcesAndSourcesCore > &post_proc_volume) const
Create finite element to post-process results.
MoFEMErrorCode getOptions()
get options
Definition: Simple.cpp:48
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:433
continuous field
Definition: definitions.h:171
MoFEMErrorCode printError(Vec ghost_vec)
Print error.
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
Definition: Simple.cpp:256
MoFEMErrorCode assembleGhostVector(Vec ghost_vec) const
Assemble error vector.
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:61
field with C-1 continuity
Definition: definitions.h:174

Variable Documentation

◆ help

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