v0.14.0
Loading...
Searching...
No Matches
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[] 
)

< 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 207 of file analytical_poisson_field_split.cpp.

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

Variable Documentation

◆ debug

const bool debug = false
static

Definition at line 21 of file analytical_poisson_field_split.cpp.

◆ help

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

Definition at line 20 of file analytical_poisson_field_split.cpp.