v0.13.0
Classes | Typedefs | Functions | Variables
scalar_check_approximation.cpp File Reference
#include <MoFEM.hpp>

Go to the source code of this file.

Classes

struct  ApproxFunctionsImpl< DIM >
 
struct  ElementsAndOps< DIM >
 
struct  ElementsAndOps< 2 >
 
struct  ElementsAndOps< 3 >
 
struct  ApproxFunctionsImpl< 2 >
 
struct  ApproxFunctionsImpl< 3 >
 
struct  OpValsDiffVals
 
struct  OpCheckValsDiffVals
 

Typedefs

using EntData = EntitiesFieldData::EntData
 
using DomainEle = ElementsAndOps< SPACE_DIM >::DomainEle
 
using DomainEleOp = ElementsAndOps< SPACE_DIM >::DomainEleOp
 
using ApproxFunctions = ApproxFunctionsImpl< SPACE_DIM >
 

Functions

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

Variables

static char help [] = "...\n\n"
 
static int approx_order = 4
 
constexpr int SPACE_DIM
 

Typedef Documentation

◆ ApproxFunctions

Definition at line 116 of file scalar_check_approximation.cpp.

◆ DomainEle

Definition at line 49 of file scalar_check_approximation.cpp.

◆ DomainEleOp

Definition at line 50 of file scalar_check_approximation.cpp.

◆ EntData

Function Documentation

◆ main()

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

Definition at line 304 of file scalar_check_approximation.cpp.

304  {
305 
306  MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
307 
308  try {
309 
310  DMType dm_name = "DMMOFEM";
311  CHKERR DMRegister_MoFEM(dm_name);
312 
313  moab::Core mb_instance;
314  moab::Interface &moab = mb_instance;
315 
316  // Add logging channel for example
317  auto core_log = logging::core::get();
318  core_log->add_sink(
319  LogManager::createSink(LogManager::getStrmWorld(), "AT"));
320  LogManager::setLog("AT");
321  MOFEM_LOG_TAG("AT", "atom_test");
322 
323  // Create MoFEM instance
324  MoFEM::Core core(moab);
325  MoFEM::Interface &m_field = core;
326 
327  Simple *simple_interface = m_field.getInterface<Simple>();
328  PipelineManager *pipeline_mng = m_field.getInterface<PipelineManager>();
329  CHKERR simple_interface->getOptions();
330  CHKERR simple_interface->loadFile("", "");
331 
332  // Declare elements
333  enum bases {
334  AINSWORTH,
335  AINSWORTH_LOBATTO,
336  DEMKOWICZ,
337  BERNSTEIN,
338  LASBASETOP
339  };
340  const char *list_bases[] = {"ainsworth", "ainsworth_labatto", "demkowicz",
341  "bernstein"};
342  PetscBool flg;
343  PetscInt choice_base_value = AINSWORTH;
344  CHKERR PetscOptionsGetEList(PETSC_NULL, NULL, "-base", list_bases,
345  LASBASETOP, &choice_base_value, &flg);
346 
347  if (flg != PETSC_TRUE)
348  SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSIBLE_CASE, "base not set");
350  if (choice_base_value == AINSWORTH)
352  if (choice_base_value == AINSWORTH_LOBATTO)
353  base = AINSWORTH_LOBATTO_BASE;
354  else if (choice_base_value == DEMKOWICZ)
355  base = DEMKOWICZ_JACOBI_BASE;
356  else if (choice_base_value == BERNSTEIN)
358 
359  enum spaces { H1SPACE, L2SPACE, LASBASETSPACE };
360  const char *list_spaces[] = {"h1", "l2"};
361  PetscInt choice_space_value = H1SPACE;
362  CHKERR PetscOptionsGetEList(PETSC_NULL, NULL, "-space", list_spaces,
363  LASBASETSPACE, &choice_space_value, &flg);
364  if (flg != PETSC_TRUE)
365  SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSIBLE_CASE, "space not set");
366  FieldSpace space = H1;
367  if (choice_space_value == H1SPACE)
368  space = H1;
369  else if (choice_space_value == L2SPACE)
370  space = L2;
371 
372  CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &approx_order,
373  PETSC_NULL);
374 
375  CHKERR simple_interface->addDomainField("FIELD1", space, base, 1);
376  CHKERR simple_interface->setFieldOrder("FIELD1", approx_order);
377 
378  Range edges, faces;
379  CHKERR moab.get_entities_by_dimension(0, 1, edges);
380  CHKERR moab.get_entities_by_dimension(0, 2, faces);
381 
382  if (choice_base_value != BERNSTEIN) {
383  Range rise_order;
384 
385  int i = 0;
386  for (auto e : edges) {
387  if (!(i % 2)) {
388  rise_order.insert(e);
389  }
390  ++i;
391  }
392 
393  for (auto f : faces) {
394  if (!(i % 3)) {
395  rise_order.insert(f);
396  }
397  ++i;
398  }
399 
400  Range rise_twice;
401  for (auto e : rise_order) {
402  if (!(i % 2)) {
403  rise_twice.insert(e);
404  }
405  ++i;
406  }
407 
408  CHKERR simple_interface->setFieldOrder("FIELD1", approx_order + 1,
409  &rise_order);
410 
411  CHKERR simple_interface->setFieldOrder("FIELD1", approx_order + 2,
412  &rise_twice);
413  }
414 
415  CHKERR simple_interface->setUp();
416  auto dm = simple_interface->getDM();
417 
418  VectorDouble vals;
419  auto jac_ptr = boost::make_shared<MatrixDouble>();
420  auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
421  auto det_ptr = boost::make_shared<VectorDouble>();
422  MatrixDouble diff_vals;
423 
424  auto assemble_matrices_and_vectors = [&]() {
426  if (SPACE_DIM == 2)
427  pipeline_mng->getOpDomainRhsPipeline().push_back(
428  new OpSetHOWeightsOnFace());
429 
432 
433  if (SPACE_DIM == 2) {
434  pipeline_mng->getOpDomainLhsPipeline().push_back(
435  new OpCalculateHOJacForFace(jac_ptr));
436  pipeline_mng->getOpDomainLhsPipeline().push_back(
437  new OpInvertMatrix<SPACE_DIM>(jac_ptr, det_ptr, inv_jac_ptr));
438  pipeline_mng->getOpDomainLhsPipeline().push_back(
439  new OpSetHOWeightsOnFace());
440  }
441 
442  if (SPACE_DIM == 3) {
443  pipeline_mng->getOpDomainLhsPipeline().push_back(
444  new OpCalculateHOJacVolume(jac_ptr));
445  pipeline_mng->getOpDomainLhsPipeline().push_back(
446  new OpInvertMatrix<SPACE_DIM>(jac_ptr, det_ptr, nullptr));
447  pipeline_mng->getOpDomainLhsPipeline().push_back(
448  new OpSetHOWeights(det_ptr));
449  pipeline_mng->getOpDomainRhsPipeline().push_back(
450  new OpCalculateHOJacVolume(jac_ptr));
451  pipeline_mng->getOpDomainRhsPipeline().push_back(
452  new OpInvertMatrix<3>(jac_ptr, det_ptr, nullptr));
453  pipeline_mng->getOpDomainRhsPipeline().push_back(
454  new OpSetHOWeights(det_ptr));
455  }
456 
459  pipeline_mng->getOpDomainLhsPipeline().push_back(new OpMass(
460  "FIELD1", "FIELD1", [](double, double, double) { return 1.; }));
461  pipeline_mng->getOpDomainRhsPipeline().push_back(
462  new OpSource("FIELD1", ApproxFunctions::fUn));
463 
464  auto integration_rule = [](int, int, int p_data) {
465  return 2 * p_data + 1;
466  };
467  CHKERR pipeline_mng->setDomainRhsIntegrationRule(integration_rule);
468  CHKERR pipeline_mng->setDomainLhsIntegrationRule(integration_rule);
469 
471  };
472 
473  auto solve_problem = [&] {
475  auto solver = pipeline_mng->createKSP();
476  CHKERR KSPSetFromOptions(solver);
477  CHKERR KSPSetUp(solver);
478 
479  auto dm = simple_interface->getDM();
480  auto D = smartCreateDMVector(dm);
481  auto F = smartVectorDuplicate(D);
482 
483  CHKERR KSPSolve(solver, F, D);
484  CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
485  CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
486  CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
487 
489  };
490 
491  auto check_solution = [&] {
493 
494  boost::shared_ptr<VectorDouble> ptr_values(new VectorDouble());
495  boost::shared_ptr<MatrixDouble> ptr_diff_vals(new MatrixDouble());
496 
497  pipeline_mng->getDomainLhsFE().reset();
498  pipeline_mng->getOpDomainRhsPipeline().clear();
499 
500  if (SPACE_DIM == 2) {
501  pipeline_mng->getOpDomainRhsPipeline().push_back(
502  new OpCalculateHOJacForFace(jac_ptr));
503  pipeline_mng->getOpDomainRhsPipeline().push_back(
504  new OpInvertMatrix<SPACE_DIM>(jac_ptr, det_ptr, inv_jac_ptr));
505  pipeline_mng->getOpDomainRhsPipeline().push_back(
506  new OpSetInvJacSpaceForFaceImpl<2>(space, inv_jac_ptr));
507  }
508 
509  if (SPACE_DIM == 3) {
510  pipeline_mng->getOpDomainRhsPipeline().push_back(
511  new OpCalculateHOJacVolume(jac_ptr));
512  pipeline_mng->getOpDomainRhsPipeline().push_back(
513  new OpInvertMatrix<SPACE_DIM>(jac_ptr, det_ptr, inv_jac_ptr));
514  pipeline_mng->getOpDomainRhsPipeline().push_back(
515  new OpSetHOInvJacToScalarBases(space, inv_jac_ptr));
516  }
517 
518  pipeline_mng->getOpDomainRhsPipeline().push_back(
519  new OpValsDiffVals(vals, diff_vals, true));
520  pipeline_mng->getOpDomainRhsPipeline().push_back(
521  new OpCalculateScalarFieldValues("FIELD1", ptr_values));
522  pipeline_mng->getOpDomainRhsPipeline().push_back(
524  ptr_diff_vals));
525  pipeline_mng->getOpDomainRhsPipeline().push_back(new OpCheckValsDiffVals(
526  vals, diff_vals, ptr_values, ptr_diff_vals, true));
527 
528  CHKERR pipeline_mng->loopFiniteElements();
529 
531  };
532 
533  CHKERR assemble_matrices_and_vectors();
534  CHKERR solve_problem();
535  CHKERR check_solution();
536  }
537  CATCH_ERRORS;
538 
540 }
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:385
FieldApproximationBase
approximation base
Definition: definitions.h:71
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:73
@ AINSWORTH_LOBATTO_BASE
Definition: definitions.h:75
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:79
@ AINSWORTH_BERNSTEIN_BEZIER_BASE
Definition: definitions.h:77
FieldSpace
approximation spaces
Definition: definitions.h:95
@ L2
field with C-1 continuity
Definition: definitions.h:101
@ H1
continuous field
Definition: definitions.h:98
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
@ MOFEM_IMPOSIBLE_CASE
Definition: definitions.h:48
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
#define CHKERR
Inline error check.
Definition: definitions.h:548
auto integration_rule
auto smartCreateDMVector
Get smart vector from DM.
Definition: DMMoFEM.hpp:956
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMMoFEM.cpp:481
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:59
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:342
FTensor::Index< 'i', SPACE_DIM > i
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:88
UBlasVector< double > VectorDouble
Definition: Types.hpp:79
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
OpCalculateHOJacForFaceImpl< 2 > OpCalculateHOJacForFace
CoreTmp< 0 > Core
Definition: Core.hpp:1096
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1965
SmartPetscObj< Vec > smartVectorDuplicate(SmartPetscObj< Vec > &vec)
Create duplicate vector of smart vector.
const double D
diffusivity
static int approx_order
static char help[]
constexpr int SPACE_DIM
static FTensor::Tensor1< double, 3 > fUn(const double x, const double y)
Core (interface) class.
Definition: Core.hpp:92
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:85
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:125
Deprecated interface functions.
Calculate jacobian on Hex or other volume which is not simplex.
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Get value at integration points for scalar field.
Set inverse jacobian to base functions.
Set inverse jacobian to base functions.
Modify integration weights on face to take in account higher-order geometry.
Transform local reference derivatives of shape functions to global derivatives.
PipelineManager interface.
Simple interface for fast problem set-up.
Definition: Simple.hpp:33
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:293
MoFEMErrorCode getOptions()
get options
Definition: Simple.cpp:216
MoFEMErrorCode getDM(DM *dm)
Get DM.
Definition: Simple.cpp:715
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name)
Load mesh file.
Definition: Simple.cpp:230
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
Definition: Simple.cpp:508
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
Definition: Simple.cpp:668
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.

Variable Documentation

◆ approx_order

int approx_order = 4
static

◆ help

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

◆ SPACE_DIM

constexpr int SPACE_DIM
constexpr
Initial value:
=
#define EXECUTABLE_DIMENSION
Definition: plastic.cpp:24
Examples
scalar_check_approximation.cpp.

Definition at line 45 of file scalar_check_approximation.cpp.