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

Go to the source code of this file.

Classes

struct  OpDivergence
 
struct  OpFlux
 

Typedefs

using FaceEle = MoFEM::FaceElementForcesAndSourcesCoreSwitch< FaceElementForcesAndSourcesCore::NO_HO_GEOMETRY|FaceElementForcesAndSourcesCore::NO_CONTRAVARIANT_TRANSFORM_HDIV|FaceElementForcesAndSourcesCore::NO_COVARIANT_TRANSFORM_HCURL >
 
using EdgeEle = MoFEM::EdgeElementForcesAndSourcesCoreSwitch< EdgeElementForcesAndSourcesCore::NO_HO_GEOMETRY|EdgeElementForcesAndSourcesCore::NO_COVARIANT_TRANSFORM_HCURL >
 
using FaceEleOp = FaceEle::UserDataOperator
 
using EdgeEleOp = EdgeEle::UserDataOperator
 

Functions

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

Variables

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

Typedef Documentation

◆ EdgeEle

using EdgeEle = MoFEM::EdgeElementForcesAndSourcesCoreSwitch< EdgeElementForcesAndSourcesCore::NO_HO_GEOMETRY | EdgeElementForcesAndSourcesCore::NO_COVARIANT_TRANSFORM_HCURL>

Definition at line 36 of file hcurl_divergence_operator_2d.cpp.

◆ EdgeEleOp

Definition at line 39 of file hcurl_divergence_operator_2d.cpp.

◆ FaceEle

using FaceEle = MoFEM::FaceElementForcesAndSourcesCoreSwitch< FaceElementForcesAndSourcesCore::NO_HO_GEOMETRY | FaceElementForcesAndSourcesCore::NO_CONTRAVARIANT_TRANSFORM_HDIV | FaceElementForcesAndSourcesCore::NO_COVARIANT_TRANSFORM_HCURL>

Definition at line 32 of file hcurl_divergence_operator_2d.cpp.

◆ FaceEleOp

using FaceEleOp = FaceEle::UserDataOperator

Function Documentation

◆ main()

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

Definition at line 96 of file hcurl_divergence_operator_2d.cpp.

96  {
97 
98  MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
99 
100  try {
101 
102  moab::Core mb_instance;
103  moab::Interface &moab = mb_instance;
104 
105  // Read mesh to MOAB
106  CHKERR moab.load_file("rectangle.h5m", 0, "");
107  ParallelComm *pcomm = ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
108  if (pcomm == NULL)
109  pcomm = new ParallelComm(&moab, PETSC_COMM_WORLD);
110 
111  // Create MoFEM instance
112  MoFEM::Core core(moab);
113  MoFEM::Interface &m_field = core;
114 
115  // set entities bit level
116  BitRefLevel bit_level0 = BitRefLevel().set(0);
117  CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevelByDim(
118  0, 2, bit_level0);
119 
120  // Declare elements
121  enum bases { AINSWORTH, DEMKOWICZ, LASBASETOP };
122  const char *list_bases[] = {"ainsworth", "demkowicz"};
123  PetscBool flg;
124  PetscInt choice_base_value = AINSWORTH;
125  CHKERR PetscOptionsGetEList(PETSC_NULL, NULL, "-base", list_bases,
126  LASBASETOP, &choice_base_value, &flg);
127  if (flg != PETSC_TRUE)
128  SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSIBLE_CASE, "base not set");
130  if (choice_base_value == AINSWORTH)
132  else if (choice_base_value == DEMKOWICZ)
133  base = DEMKOWICZ_JACOBI_BASE;
134  int order = 1;
135  CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &order, PETSC_NULL);
136 
137  CHKERR m_field.add_field("FIELD1", HCURL, base, 1);
138  CHKERR m_field.add_finite_element("FACE_FE");
139  // Define rows/cols and element data
140  CHKERR m_field.modify_finite_element_add_field_row("FACE_FE", "FIELD1");
141  CHKERR m_field.modify_finite_element_add_field_col("FACE_FE", "FIELD1");
142  CHKERR m_field.modify_finite_element_add_field_data("FACE_FE", "FIELD1");
143 
144  CHKERR m_field.add_finite_element("EDGE_FE");
145  // Define rows/cols and element data
146  CHKERR m_field.modify_finite_element_add_field_row("EDGE_FE", "FIELD1");
147  CHKERR m_field.modify_finite_element_add_field_col("EDGE_FE", "FIELD1");
148  CHKERR m_field.modify_finite_element_add_field_data("EDGE_FE", "FIELD1");
149 
150  // Problem
151  CHKERR m_field.add_problem("TEST_PROBLEM");
152  // set finite elements for problem
153  CHKERR m_field.modify_problem_add_finite_element("TEST_PROBLEM", "FACE_FE");
154  CHKERR m_field.modify_problem_add_finite_element("TEST_PROBLEM", "EDGE_FE");
155  // set refinement level for problem
156  CHKERR m_field.modify_problem_ref_level_add_bit("TEST_PROBLEM", bit_level0);
157 
158  // Add entities
159 
160  CHKERR m_field.add_ents_to_field_by_type(0, MBTRI, "FIELD1");
161  // Set order
162  CHKERR m_field.set_field_order(0, MBTRI, "FIELD1", order);
163  CHKERR m_field.set_field_order(0, MBEDGE, "FIELD1", order);
164 
165  // Add entities to elements
166  CHKERR m_field.add_ents_to_finite_element_by_type(0, MBTRI, "FACE_FE");
167 
168  auto set_edge_elements_entities_on_mesh_skin = [&]() {
170  Range faces;
171  CHKERR moab.get_entities_by_type(0, MBTRI, faces, false);
172  Skinner skin(&m_field.get_moab());
173  Range faces_skin;
174  CHKERR skin.find_skin(0, faces, false, faces_skin);
175  CHKERR m_field.add_ents_to_finite_element_by_type(faces_skin, MBEDGE,
176  "EDGE_FE");
178  };
179  CHKERR set_edge_elements_entities_on_mesh_skin();
180 
181  // Build database
182  CHKERR m_field.build_fields();
183  // build finite elemnts
184  CHKERR m_field.build_finite_elements();
185  // build adjacencies
186  CHKERR m_field.build_adjacencies(bit_level0);
187 
188  // build problem
189  ProblemsManager *prb_mng_ptr;
190  CHKERR m_field.getInterface(prb_mng_ptr);
191  CHKERR prb_mng_ptr->buildProblem("TEST_PROBLEM", true);
192  // Partition
193  CHKERR prb_mng_ptr->partitionSimpleProblem("TEST_PROBLEM");
194  CHKERR prb_mng_ptr->partitionFiniteElements("TEST_PROBLEM");
195  CHKERR prb_mng_ptr->partitionGhostDofs("TEST_PROBLEM");
196 
197  // integration rule
198  auto rule = [&](int, int, int p) { return p; };
199 
200  auto calculate_divergence = [&]() {
201  double div = 0;
202  FaceEle fe_face(m_field);
203  fe_face.getRuleHook = rule;
204  MatrixDouble inv_jac(2, 2), jac(2, 2);
205  fe_face.getOpPtrVector().push_back(new OpCalculateJacForFace(jac));
206  fe_face.getOpPtrVector().push_back(new OpCalculateInvJacForFace(inv_jac));
207  fe_face.getOpPtrVector().push_back(new OpMakeHdivFromHcurl());
208  fe_face.getOpPtrVector().push_back(
210  fe_face.getOpPtrVector().push_back(new OpSetInvJacHcurlFace(inv_jac));
211  fe_face.getOpPtrVector().push_back(new OpDivergence(div));
212  CHKERR m_field.loop_finite_elements("TEST_PROBLEM", "FACE_FE", fe_face);
213  return div;
214  };
215 
216  auto calculate_flux = [&]() {
217  double flux = 0;
218  EdgeEle fe_edge(m_field);
219  fe_edge.getRuleHook = rule;
220  fe_edge.getOpPtrVector().push_back(
222  fe_edge.getOpPtrVector().push_back(new OpFlux(flux));
223  CHKERR m_field.loop_finite_elements("TEST_PROBLEM", "EDGE_FE", fe_edge);
224  return flux;
225  };
226 
227  const double div = calculate_divergence();
228  const double flux = calculate_flux();
229 
230  PetscPrintf(PETSC_COMM_WORLD, "Div = %4.3e Flux = %3.4e Error = %4.3e\n",
231  div, flux, div + flux);
232 
233  constexpr double tol = 1e-8;
234  if (std::abs(div + flux) > tol)
235  SETERRQ2(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
236  "Test failed (div != flux) %3.4e != %3.4e", div, flux);
237  }
238  CATCH_ERRORS;
239 
241 }
virtual MoFEMErrorCode modify_finite_element_add_field_col(const std::string &fe_name, const std::string &name_row)=0
set field col which finite element use
virtual MoFEMErrorCode add_ents_to_finite_element_by_type(const EntityHandle entities, const EntityType type, const std::string &name, const bool recursive=true)=0
add entities to finite element
Deprecated interface functions.
Problem manager is used to build and partition problems.
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
virtual moab::Interface & get_moab()=0
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:74
virtual MoFEMErrorCode modify_problem_add_finite_element(const std::string &name_problem, const std::string &fe_name)=0
add finite element to problem, this add entities assigned to finite element to a particular problem
MoFEMErrorCode buildProblem(const std::string name, const bool square_matrix, int verb=VERBOSE)
build problem data structures
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:482
virtual MoFEMErrorCode add_field(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_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add field.
virtual MoFEMErrorCode set_field_order(const EntityHandle meshset, const EntityType type, const std::string &name, const ApproximationOrder order, int verb=DEFAULT_VERBOSITY)=0
Set order approximation of the entities in the field.
Core (interface) class.
Definition: Core.hpp:50
Edge finite elementUser is implementing own operator at Gauss points level, by own object derived fro...
static char help[]
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 getInterface(const MOFEMuuid &uuid, IFACE *&iface) const
Get interface by uuid and return reference to pointer of interface.
Calculate inverse of jacobian for face element.
double tol
MoFEM::FaceElementForcesAndSourcesCore FaceEle
FieldApproximationBase
approximation base
Definition: definitions.h:149
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
constexpr int order
virtual MoFEMErrorCode add_ents_to_field_by_type(const Range &ents, const EntityType type, const std::string &name, int verb=DEFAULT_VERBOSITY)=0
Add entities to field meshset.
Managing BitRefLevels.
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:151
virtual MoFEMErrorCode modify_problem_ref_level_add_bit(const std::string &name_problem, const BitRefLevel &bit)=0
add ref level to problem
field with continuous tangents
Definition: definitions.h:177
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
#define CHKERR
Inline error check.
Definition: definitions.h:601
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string &name_filed)=0
set finite element field data
MoFEMErrorCode partitionGhostDofs(const std::string name, int verb=VERBOSE)
determine ghost nodes
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:290
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:51
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
brief Transform local reference derivatives of shape function to global derivatives for face
virtual MoFEMErrorCode add_problem(const std::string &name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add problem.
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1791
MoFEMErrorCode partitionFiniteElements(const std::string name, bool part_from_moab=false, int low_proc=-1, int hi_proc=-1, int verb=VERBOSE)
partition finite elementsFunction which partition finite elements based on dofs partitioning....
Apply contravariant (Piola) transfer to Hdiv space on face.
virtual MoFEMErrorCode modify_finite_element_add_field_row(const std::string &fe_name, const std::string &name_row)=0
set field row which finite element use
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:412
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:438
Calculate jacobian for face element.
MoFEMErrorCode partitionSimpleProblem(const std::string name, int verb=VERBOSE)
partition problem dofs
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elementsBuild finite element data structures. Have to be run before problem and adjacenc...
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:61
DEPRECATED MoFEMErrorCode loop_finite_elements(const Problem *problem_ptr, const std::string &fe_name, FEMethod &method, int lower_rank, int upper_rank, MoFEMTypes bh, int verb=DEFAULT_VERBOSITY)
Make Hdiv space from Hcurl space in 2d.

Variable Documentation

◆ help

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