v0.9.0
Classes | Functions | Variables
analytical_poisson_field_split.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  OpS
 

Functions

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

Variables

static char help [] = "...\n\n"
 
static const bool debug = false
 

Function Documentation

◆ main()

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

Definition at line 219 of file analytical_poisson_field_split.cpp.

219  {
220 
221  // Initialize PETSc
222  MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
223  // Create MoAB database
224  moab::Core moab_core; // create database
225  moab::Interface &moab = moab_core; // create interface to database
226 
227  try {
228 
229  // Get command line options
230  int order = 3; // default approximation order
231  PetscBool flg_test = PETSC_FALSE; // true check if error is numerical error
232  CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "Poisson's problem options",
233  "none");
234  // Set approximation order
235  CHKERR PetscOptionsInt("-order", "approximation order", "", order, &order,
236  PETSC_NULL);
237  // Set testing (used by CTest)
238  CHKERR PetscOptionsBool("-test", "if true is ctest", "", flg_test,
239  &flg_test, PETSC_NULL);
240  ierr = PetscOptionsEnd();
241  CHKERRG(ierr);
242 
243  // Create MoFEM database and link it to MoAB
244  MoFEM::Core mofem_core(moab); // create database
245  MoFEM::Interface &m_field = mofem_core; // create interface to database
246  // Register DM Manager
247  CHKERR DMRegister_MoFEM("DMMOFEM"); // register MoFEM DM in PETSc
248 
249  // Create vector to store approximation global error
250  Vec global_error;
251  CHKERR PoissonExample::AuxFunctions(m_field).createGhostVec(&global_error);
252 
253  // First we crate elements, implementation of elements is problem
254  // independent, we do not know yet what fields are present in the problem,
255  // or even we do not decided yet what approximation base or spaces we are
256  // going to use. Implementation of element is free from those constrains and
257  // can be used in different context.
258 
259  // Elements used by KSP & DM to assemble system of equations
260  boost::shared_ptr<ForcesAndSourcesCore>
261  domain_lhs_fe; ///< Volume element for the matrix
262  boost::shared_ptr<ForcesAndSourcesCore>
263  boundary_lhs_fe; ///< Boundary element for the matrix
264  boost::shared_ptr<ForcesAndSourcesCore>
265  domain_rhs_fe; ///< Volume element to assemble vector
266  boost::shared_ptr<ForcesAndSourcesCore>
267  boundary_rhs_fe; ///< Volume element to assemble vector
268  boost::shared_ptr<ForcesAndSourcesCore>
269  domain_error; ///< Volume element evaluate error
270  boost::shared_ptr<ForcesAndSourcesCore>
271  post_proc_volume; ///< Volume element to Post-process results
272  boost::shared_ptr<ForcesAndSourcesCore> null; ///< Null element do nothing
273  boost::shared_ptr<ForcesAndSourcesCore> boundary_penalty_lhs_fe;
274  {
275  // Add problem specific operators the generic finite elements to calculate
276  // matrices and vectors.
279  ExactFunction(), ExactLaplacianFunction(), domain_lhs_fe,
280  boundary_lhs_fe, domain_rhs_fe, boundary_rhs_fe, false);
281  // Add problem specific operators the generic finite elements to calculate
282  // error on elements and global error in H1 norm
285  global_error, domain_error);
286  // Post-process results
288  .creatFEToPostProcessResults(post_proc_volume);
289 
290  const double beta = 1;
291  boundary_penalty_lhs_fe = boost::shared_ptr<ForcesAndSourcesCore>(
292  new FaceElementForcesAndSourcesCore(m_field));
293  boundary_penalty_lhs_fe->getRuleHook = PoissonExample::FaceRule();
294  boundary_penalty_lhs_fe->getOpPtrVector().push_back(new OpS(beta));
295  boundary_rhs_fe->getOpPtrVector().push_back(
296  new PoissonExample::Op_g(ExactFunction(), "U", beta));
297  }
298 
299  // Get simple interface is simplified version enabling quick and
300  // easy construction of problem.
301  Simple *simple_interface;
302  // Query interface and get pointer to Simple interface
303  CHKERR m_field.getInterface(simple_interface);
304 
305  // Build problem with simple interface
306  {
307 
308  // Get options for simple interface from command line
309  CHKERR simple_interface->getOptions();
310  // Load mesh file to database
311  CHKERR simple_interface->loadFile();
312 
313  // Add field on domain and boundary. Field is declared by space and base
314  // and rank. space can be H1. Hcurl, Hdiv and L2, base can be
315  // AINSWORTH_LEGENDRE_BASE, DEMKOWICZ_JACOBI_BASE and more, where rank is
316  // number of coefficients for dof.
317  //
318  // Simple interface assumes that there are four types of field; 1) defined
319  // on body domain, 2) fields defined on body boundary, 3) skeleton field
320  // defined on finite element skeleton. Finally data field is defined on
321  // body domain. Data field is not solved but used for post-process or to
322  // keep material parameters, etc. Here we using it to calculate
323  // approximation error on elements.
324 
325  // Add domain filed "U" in space H1 and Legendre base, Ainsworth recipe is
326  // used to construct base functions.
327  CHKERR simple_interface->addDomainField("U", H1, AINSWORTH_LEGENDRE_BASE,
328  1);
329  // Add Lagrange multiplier field on body boundary
330  CHKERR simple_interface->addBoundaryField("L", H1,
332  // Add error (data) field, we need only L2 norm. Later order is set to 0,
333  // so this is piecewise discontinuous constant approx., i.e. 1 DOF for
334  // element. You can use more DOFs and collate moments of error to drive
335  // anisotropic h/p-adaptivity, however this is beyond this example.
336  CHKERR simple_interface->addDataField("ERROR", L2,
338 
339  // Set fields order domain and boundary fields.
340  CHKERR simple_interface->setFieldOrder("U",
341  order); // to approximate function
342  CHKERR simple_interface->setFieldOrder("L",
343  order); // to Lagrange multipliers
344  CHKERR simple_interface->setFieldOrder(
345  "ERROR", 0); // approximation order for error
346 
347  // Setup problem. At that point database is constructed, i.e. fields,
348  // finite elements and problem data structures with local and global
349  // indexing.
350  CHKERR simple_interface->setUp();
351  }
352 
353  // Get access to PETSC-MoFEM DM manager.
354  // or more derails see
355  // <http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/DM/index.html>
356  // Form that point internal MoFEM data structures are linked with PETSc
357  // interface. If DM functions contains string MoFEM is is MoFEM specific DM
358  // interface function, otherwise DM function part of standard PETSc
359  // interface.
360 
361  DM dm;
362  // Get dm
363  CHKERR simple_interface->getDM(&dm);
364 
365  // Solve problem, only PETEc interface here.
366  {
367 
368  // Create the right hand side vector and vector of unknowns
369  Vec F, D;
370  CHKERR DMCreateGlobalVector(dm, &F);
371  // Create unknown vector by creating duplicate copy of F vector. only
372  // structure is duplicated no values.
373  CHKERR VecDuplicate(F, &D);
374 
375  DM dm_sub_KK, dm_sub_LU;
376  ublas::matrix<Mat> nested_matrices(2, 2);
377  ublas::vector<IS> nested_is(2);
378 
379  CHKERR DMCreate(PETSC_COMM_WORLD, &dm_sub_KK);
380  CHKERR DMSetType(dm_sub_KK, "DMMOFEM");
381  CHKERR DMMoFEMCreateSubDM(dm_sub_KK, dm, "SUB_KK");
382  CHKERR DMSetFromOptions(dm_sub_KK);
383  CHKERR DMMoFEMSetSquareProblem(dm_sub_KK, PETSC_TRUE);
384  CHKERR DMMoFEMAddSubFieldRow(dm_sub_KK, "U");
385  CHKERR DMMoFEMAddSubFieldCol(dm_sub_KK, "U");
386  CHKERR DMMoFEMAddElement(dm_sub_KK,
387  simple_interface->getDomainFEName().c_str());
388  CHKERR DMMoFEMAddElement(dm_sub_KK,
389  simple_interface->getBoundaryFEName().c_str());
390  CHKERR DMSetUp(dm_sub_KK);
391  CHKERR DMMoFEMGetSubRowIS(dm_sub_KK, &nested_is[0]);
392  CHKERR DMCreateMatrix(dm_sub_KK, &nested_matrices(0, 0));
393  domain_lhs_fe->ksp_B = nested_matrices(0, 0);
395  dm_sub_KK, simple_interface->getDomainFEName(), domain_lhs_fe);
396  boundary_penalty_lhs_fe->ksp_B = nested_matrices(0, 0);
398  simple_interface->getBoundaryFEName(),
399  boundary_penalty_lhs_fe);
400  CHKERR MatAssemblyBegin(nested_matrices(0, 0), MAT_FINAL_ASSEMBLY);
401  CHKERR MatAssemblyEnd(nested_matrices(0, 0), MAT_FINAL_ASSEMBLY);
402  CHKERR DMDestroy(&dm_sub_KK);
403  //
404  CHKERR DMCreate(PETSC_COMM_WORLD, &dm_sub_LU);
405  CHKERR DMSetType(dm_sub_LU, "DMMOFEM");
406  CHKERR DMMoFEMCreateSubDM(dm_sub_LU, dm, "SUB_LU");
407  CHKERR DMSetFromOptions(dm_sub_LU);
408  CHKERR DMMoFEMSetSquareProblem(dm_sub_LU, PETSC_FALSE);
409  CHKERR DMMoFEMAddSubFieldRow(dm_sub_LU, "L");
410  CHKERR DMMoFEMAddSubFieldCol(dm_sub_LU, "U");
411  CHKERR DMMoFEMAddElement(dm_sub_LU,
412  simple_interface->getBoundaryFEName().c_str());
413  CHKERR DMSetUp(dm_sub_LU);
414  CHKERR DMMoFEMGetSubRowIS(dm_sub_LU, &nested_is[1]);
415  CHKERR DMCreateMatrix(dm_sub_LU, &nested_matrices(1, 0));
416  boundary_lhs_fe->ksp_B = nested_matrices(1, 0);
418  dm_sub_LU, simple_interface->getBoundaryFEName(), boundary_lhs_fe);
419  CHKERR MatAssemblyBegin(nested_matrices(1, 0), MAT_FINAL_ASSEMBLY);
420  CHKERR MatAssemblyEnd(nested_matrices(1, 0), MAT_FINAL_ASSEMBLY);
421  // CHKERR MatCreateTranspose(nested_matrices(1,0),&nested_matrices(0,1));
422  CHKERR MatTranspose(nested_matrices(1, 0), MAT_INITIAL_MATRIX,
423  &nested_matrices(0, 1));
424  CHKERR DMDestroy(&dm_sub_LU);
425 
426  domain_rhs_fe->ksp_f = F;
427  CHKERR DMoFEMLoopFiniteElements(dm, simple_interface->getDomainFEName(),
428  domain_rhs_fe);
429  boundary_rhs_fe->ksp_f = F;
430  CHKERR DMoFEMLoopFiniteElements(dm, simple_interface->getBoundaryFEName(),
431  boundary_rhs_fe);
432  CHKERR VecAssemblyBegin(F);
433  CHKERR VecAssemblyEnd(F);
434 
435  Mat A;
436  nested_matrices(1, 1) = PETSC_NULL;
437 
438  if (debug) {
439  MatType type;
440  MatGetType(nested_matrices(0, 0), &type);
441  cerr << "K " << type << endl;
442  MatGetType(nested_matrices(1, 0), &type);
443  cerr << "C " << type << endl;
444  MatGetType(nested_matrices(0, 1), &type);
445  cerr << "CT " << type << endl;
446  std::string wait;
447  cerr << "UU" << endl;
448  MatView(nested_matrices(0, 0), PETSC_VIEWER_DRAW_WORLD);
449  std::cin >> wait;
450  cerr << "LU" << endl;
451  MatView(nested_matrices(1, 0), PETSC_VIEWER_DRAW_WORLD);
452  std::cin >> wait;
453  cerr << "UL" << endl;
454  MatView(nested_matrices(0, 1), PETSC_VIEWER_DRAW_WORLD);
455  std::cin >> wait;
456  }
457 
458  CHKERR MatCreateNest(PETSC_COMM_WORLD, 2, &nested_is[0], 2, &nested_is[0],
459  &nested_matrices(0, 0), &A);
460 
461  // Create solver and link it to DM
462  KSP solver;
463  CHKERR KSPCreate(PETSC_COMM_WORLD, &solver);
464  CHKERR KSPSetFromOptions(solver);
465  // Set operators
466  CHKERR KSPSetOperators(solver, A, A);
467  PC pc;
468  CHKERR KSPGetPC(solver, &pc);
469  PetscBool is_pcfs = PETSC_FALSE;
470  PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
471  if (is_pcfs) {
472  CHKERR PCFieldSplitSetIS(pc, NULL, nested_is[0]);
473  CHKERR PCFieldSplitSetIS(pc, NULL, nested_is[1]);
474  } else {
475  SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
476  "Only works with pre-conditioner PCFIELDSPLIT");
477  }
478  // Set-up solver, is type of solver and pre-conditioners
479  CHKERR KSPSetUp(solver);
480  // At solution process, KSP solver using DM creates matrices, Calculate
481  // values of the left hand side and the right hand side vector. then
482  // solves system of equations. Results are stored in vector D.
483  CHKERR KSPSolve(solver, F, D);
484 
485  // Scatter solution on the mesh. Stores unknown vector on field on the
486  // mesh.
487  CHKERR DMoFEMMeshToGlobalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
488 
489  // Clean data. Solver and vector are not needed any more.
490  CHKERR KSPDestroy(&solver);
491  for (int i = 0; i != 2; i++) {
492  CHKERR ISDestroy(&nested_is[i]);
493  for (int j = 0; j != 2; j++) {
494  CHKERR MatDestroy(&nested_matrices(i, j));
495  }
496  }
497  CHKERR MatDestroy(&A);
498  CHKERR VecDestroy(&D);
499  CHKERR VecDestroy(&F);
500  }
501 
502  // Calculate error
503  {
504  // Loop over all elements in mesh, and run users operators on each
505  // element.
506  CHKERR DMoFEMLoopFiniteElements(dm, simple_interface->getDomainFEName(),
507  domain_error);
509  global_error);
510  CHKERR PoissonExample::AuxFunctions(m_field).printError(global_error);
511  if (flg_test == PETSC_TRUE) {
512  CHKERR PoissonExample::AuxFunctions(m_field).testError(global_error);
513  }
514  }
515 
516  {
517  // Loop over all elements in the mesh and for each execute
518  // post_proc_volume element and operators on it.
519  CHKERR DMoFEMLoopFiniteElements(dm, simple_interface->getDomainFEName(),
520  post_proc_volume);
521  // Write results
522  CHKERR boost::static_pointer_cast<PostProcVolumeOnRefinedMesh>(
523  post_proc_volume)
524  ->writeFile("out_vol.h5m");
525  }
526 
527  // Destroy DM, no longer needed.
528  CHKERR DMDestroy(&dm);
529 
530  // Destroy ghost vector
531  CHKERR VecDestroy(&global_error);
532  }
533  CATCH_ERRORS;
534 
535  // finish work cleaning memory, getting statistics, etc.
537 
538  return 0;
539 }
Assemble constrains vector.
PetscErrorCode DMMoFEMGetSubRowIS(DM dm, IS *is)
get sub problem is
Definition: DMMMoFEM.cpp:242
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
Definition: Simple.cpp:507
const std::string getBoundaryFEName() const
Definition: Simple.hpp:222
MoFEMErrorCode getDM(DM *dm)
Get DM.
Definition: Simple.cpp:517
PetscErrorCode DMMoFEMSetSquareProblem(DM dm, PetscBool square_problem)
set squared problemIt if true is assumed that matrix has the same indexing on rows and columns....
Definition: DMMMoFEM.cpp:367
PetscErrorCode DMMoFEMAddSubFieldCol(DM dm, const char field_name[], EntityType lo_type=MBVERTEX, EntityType hi_type=MBMAXTYPE)
Definition: DMMMoFEM.cpp:210
PetscErrorCode DMMoFEMAddElement(DM dm, const char fe_name[])
add element to dm
Definition: DMMMoFEM.cpp:414
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
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method)
Executes FEMethod for finite elements in DM.
Definition: DMMMoFEM.cpp:496
static const bool debug
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
PetscErrorCode DMMoFEMCreateSubDM(DM subdm, DM dm, const char problem_name[])
Must be called by user to set Sub DM MoFEM data structures.
Definition: DMMMoFEM.cpp:166
PetscErrorCode DMMoFEMAddSubFieldRow(DM dm, const char field_name[], EntityType lo_type=MBVERTEX, EntityType hi_type=MBMAXTYPE)
Definition: DMMMoFEM.cpp:187
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.
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
Set integration rule to boundary elements.
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.
#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
FaceElementForcesAndSourcesCoreSwitch< 0 > FaceElementForcesAndSourcesCore
Face finite element default.
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 char help[]
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

◆ debug

const bool debug = false
static

◆ help

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