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

Go to the source code of this file.

Classes

struct  MyOp
 
struct  CallingOp
 
struct  MyOp2
 

Typedefs

typedef tee_device< std::ostream, std::ofstream > TeeDevice
 
typedef stream< TeeDeviceTeeStream
 

Functions

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

Variables

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

Typedef Documentation

◆ TeeDevice

typedef tee_device<std::ostream, std::ofstream> TeeDevice

◆ TeeStream

typedef stream<TeeDevice> TeeStream

Function Documentation

◆ main()

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

Definition at line 180 of file forces_and_sources_testing_contact_prism_element.cpp.

180  {
181 
182  MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
183 
184  try {
185 
186  moab::Core mb_instance;
187  moab::Interface &moab = mb_instance;
188  ParallelComm *pcomm = ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
189  if (pcomm == NULL)
190  pcomm = new ParallelComm(&moab, PETSC_COMM_WORLD);
191 
192  PetscBool flg = PETSC_TRUE;
193  char mesh_file_name[255];
194 #if PETSC_VERSION_GE(3, 6, 4)
195  CHKERR PetscOptionsGetString(PETSC_NULL, "", "-my_file", mesh_file_name,
196  255, &flg);
197 #else
198  CHKERR PetscOptionsGetString(PETSC_NULL, PETSC_NULL, "-my_file",
199  mesh_file_name, 255, &flg);
200 #endif
201  if (flg != PETSC_TRUE)
202  SETERRQ(PETSC_COMM_SELF, 1, "error -my_file (mesh file not given");
203 
204  const char *option;
205  option = "";
206  CHKERR moab.load_file(mesh_file_name, 0, option);
207 
208  // Create MoFEM (Joseph) database
209  MoFEM::Core core(moab);
210  MoFEM::Interface &m_field = core;
211  PrismInterface *interface;
212  CHKERR m_field.getInterface(interface);
213 
214  // set entities bit level
215  CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevelByDim(
216  0, 3, BitRefLevel().set(0));
217  std::vector<BitRefLevel> bit_levels;
218  bit_levels.push_back(BitRefLevel().set(0));
219 
220  int ll = 1;
221  for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field, SIDESET, cit)) {
222 
223  CHKERR PetscPrintf(PETSC_COMM_WORLD, "Insert Interface %d\n",
224  cit->getMeshsetId());
225  EntityHandle cubit_meshset = cit->getMeshset();
226  {
227  // get tet enties form back bit_level
228  EntityHandle ref_level_meshset = 0;
229  CHKERR moab.create_meshset(MESHSET_SET, ref_level_meshset);
231  ->getEntitiesByTypeAndRefLevel(bit_levels.back(),
232  BitRefLevel().set(), MBTET,
233  ref_level_meshset);
235  ->getEntitiesByTypeAndRefLevel(bit_levels.back(),
236  BitRefLevel().set(), MBPRISM,
237  ref_level_meshset);
238  Range ref_level_tets;
239  CHKERR moab.get_entities_by_handle(ref_level_meshset, ref_level_tets,
240  true);
241  // get faces and test to split
242  CHKERR interface->getSides(cubit_meshset, bit_levels.back(), true, 0);
243  // set new bit level
244  bit_levels.push_back(BitRefLevel().set(ll++));
245  // split faces and
246  CHKERR interface->splitSides(ref_level_meshset, bit_levels.back(),
247  cubit_meshset, true, true, 0);
248  // clean meshsets
249  CHKERR moab.delete_entities(&ref_level_meshset, 1);
250  }
251  // update cubit meshsets
252  for (_IT_CUBITMESHSETS_FOR_LOOP_(m_field, ciit)) {
253  EntityHandle cubit_meshset = ciit->meshset;
255  ->updateMeshsetByEntitiesChildren(cubit_meshset, bit_levels.back(),
256  cubit_meshset, MBVERTEX, true);
258  ->updateMeshsetByEntitiesChildren(cubit_meshset, bit_levels.back(),
259  cubit_meshset, MBEDGE, true);
261  ->updateMeshsetByEntitiesChildren(cubit_meshset, bit_levels.back(),
262  cubit_meshset, MBTRI, true);
264  ->updateMeshsetByEntitiesChildren(cubit_meshset, bit_levels.back(),
265  cubit_meshset, MBTET, true);
266  }
267  }
268 
269  // Fields
270  CHKERR m_field.add_field("FIELD1", H1, AINSWORTH_LEGENDRE_BASE, 3);
271  CHKERR m_field.add_field("MESH_NODE_POSITIONS", H1, AINSWORTH_LEGENDRE_BASE,
272  3);
273  CHKERR m_field.add_field("FIELD2", NOFIELD, NOBASE, 3);
274 
275  auto set_no_field_vertex = [&]() {
277  // Creating and adding no field entities.
278  const double coords[] = {0, 0, 0};
279  EntityHandle no_field_vertex;
280  CHKERR m_field.get_moab().create_vertex(coords, no_field_vertex);
281  Range range_no_field_vertex;
282  range_no_field_vertex.insert(no_field_vertex);
283  CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevel(
284  range_no_field_vertex, BitRefLevel().set());
285  EntityHandle meshset = m_field.get_field_meshset("FIELD2");
286  CHKERR m_field.get_moab().add_entities(meshset, range_no_field_vertex);
288  };
289 
290  CHKERR set_no_field_vertex();
291 
292  // FE
293  CHKERR m_field.add_finite_element("TEST_FE1");
294  CHKERR m_field.add_finite_element("TEST_FE2");
295 
296  // Define rows/cols and element data
297  CHKERR m_field.modify_finite_element_add_field_row("TEST_FE1", "FIELD1");
298  CHKERR m_field.modify_finite_element_add_field_col("TEST_FE1", "FIELD1");
299  CHKERR m_field.modify_finite_element_add_field_data("TEST_FE1", "FIELD1");
300  CHKERR m_field.modify_finite_element_add_field_data("TEST_FE1",
301  "MESH_NODE_POSITIONS");
302 
303  CHKERR m_field.modify_finite_element_add_field_row("TEST_FE2", "FIELD1");
304  CHKERR m_field.modify_finite_element_add_field_col("TEST_FE2", "FIELD2");
305  CHKERR m_field.modify_finite_element_add_field_data("TEST_FE2", "FIELD1");
306  CHKERR m_field.modify_finite_element_add_field_data("TEST_FE2", "FIELD2");
307 
308  // Problem
309  CHKERR m_field.add_problem("TEST_PROBLEM");
310 
311  // set finite elements for problem
312  CHKERR m_field.modify_problem_add_finite_element("TEST_PROBLEM",
313  "TEST_FE1");
314  CHKERR m_field.modify_problem_add_finite_element("TEST_PROBLEM",
315  "TEST_FE2");
316  // set refinement level for problem
317  CHKERR m_field.modify_problem_ref_level_add_bit("TEST_PROBLEM",
318  bit_levels.back());
319 
320  // meshset consisting all entities in mesh
321  EntityHandle root_set = moab.get_root_set();
322  // add entities to field
323  CHKERR m_field.add_ents_to_field_by_type(root_set, MBTET, "FIELD1");
324  CHKERR m_field.add_ents_to_field_by_type(root_set, MBTET,
325  "MESH_NODE_POSITIONS");
326  // add entities to finite element
327  CHKERR m_field.add_ents_to_finite_element_by_type(root_set, MBPRISM,
328  "TEST_FE1", 10);
329  CHKERR m_field.add_ents_to_finite_element_by_type(root_set, MBPRISM,
330  "TEST_FE2", 10);
331 
332  // set app. order
333  // see Hierarchic Finite Element Bases on Unstructured Tetrahedral Meshes
334  // (Mark Ainsworth & Joe Coyle)
335  int order = 3;
336  CHKERR m_field.set_field_order(root_set, MBTET, "FIELD1", order);
337  CHKERR m_field.set_field_order(root_set, MBTRI, "FIELD1", order);
338  CHKERR m_field.set_field_order(root_set, MBEDGE, "FIELD1", order);
339  CHKERR m_field.set_field_order(root_set, MBVERTEX, "FIELD1", 1);
340 
341  CHKERR m_field.set_field_order(root_set, MBTET, "MESH_NODE_POSITIONS", 2);
342  CHKERR m_field.set_field_order(root_set, MBTRI, "MESH_NODE_POSITIONS", 2);
343  CHKERR m_field.set_field_order(root_set, MBEDGE, "MESH_NODE_POSITIONS", 2);
344  CHKERR m_field.set_field_order(root_set, MBVERTEX, "MESH_NODE_POSITIONS",
345  1);
346 
347  /****/
348  // build database
349  // build field
350  CHKERR m_field.build_fields();
351  // set FIELD1 from positions of 10 node tets
352  Projection10NodeCoordsOnField ent_method_field1(m_field, "FIELD1");
353  CHKERR m_field.loop_dofs("FIELD1", ent_method_field1);
354  Projection10NodeCoordsOnField ent_method_mesh_positions(
355  m_field, "MESH_NODE_POSITIONS");
356  CHKERR m_field.loop_dofs("MESH_NODE_POSITIONS", ent_method_mesh_positions);
357 
358  // build finite elemnts
359  CHKERR m_field.build_finite_elements();
360  // build adjacencies
361  CHKERR m_field.build_adjacencies(bit_levels.back());
362  // build problem
363  ProblemsManager *prb_mng_ptr;
364  CHKERR m_field.getInterface(prb_mng_ptr);
365  CHKERR prb_mng_ptr->buildProblem("TEST_PROBLEM", false);
366 
367  /****/
368  // mesh partitioning
369  // partition
370  CHKERR prb_mng_ptr->partitionSimpleProblem("TEST_PROBLEM");
371  CHKERR prb_mng_ptr->partitionFiniteElements("TEST_PROBLEM");
372  // what are ghost nodes, see Petsc Manual
373  CHKERR prb_mng_ptr->partitionGhostDofs("TEST_PROBLEM");
374 
375 
376  std::ofstream ofs("forces_and_sources_testing_contact_prism_element.txt");
377  TeeDevice my_tee(std::cout, ofs);
378  TeeStream my_split(my_tee);
379 
381  using ContactDataOp =
383 
385  fe1.getOpPtrVector().push_back(
386  new MyOp(my_split, UMDataOp::OPROW, ContactDataOp::FACEMASTER));
387  fe1.getOpPtrVector().push_back(
388  new MyOp(my_split, UMDataOp::OPROW, ContactDataOp::FACESLAVE));
389  fe1.getOpPtrVector().push_back(new MyOp(my_split, UMDataOp::OPROWCOL,
390  ContactDataOp::FACEMASTERMASTER));
391  fe1.getOpPtrVector().push_back(
392  new MyOp(my_split, UMDataOp::OPROWCOL, ContactDataOp::FACEMASTERSLAVE));
393  fe1.getOpPtrVector().push_back(
394  new MyOp(my_split, UMDataOp::OPROWCOL, ContactDataOp::FACESLAVEMASTER));
395  fe1.getOpPtrVector().push_back(
396  new MyOp(my_split, UMDataOp::OPROWCOL, ContactDataOp::FACESLAVESLAVE));
397  fe1.getOpPtrVector().push_back(new CallingOp(my_split, UMDataOp::OPCOL));
398  fe1.getOpPtrVector().push_back(new CallingOp(my_split, UMDataOp::OPROW));
399  fe1.getOpPtrVector().push_back(new CallingOp(my_split, UMDataOp::OPROWCOL));
400  CHKERR m_field.loop_finite_elements("TEST_PROBLEM", "TEST_FE1", fe1);
401 
403  fe2.getOpPtrVector().push_back(
404  new MyOp2(my_split, UMDataOp::OPCOL, ContactDataOp::FACEMASTER));
405  fe2.getOpPtrVector().push_back(
406  new MyOp2(my_split, UMDataOp::OPCOL, ContactDataOp::FACESLAVE));
407  fe2.getOpPtrVector().push_back(new MyOp2(my_split, UMDataOp::OPROWCOL,
408  ContactDataOp::FACEMASTERMASTER));
409  fe2.getOpPtrVector().push_back(new MyOp2(my_split, UMDataOp::OPROWCOL,
410  ContactDataOp::FACEMASTERSLAVE));
411  fe2.getOpPtrVector().push_back(new MyOp2(my_split, UMDataOp::OPROWCOL,
412  ContactDataOp::FACESLAVEMASTER));
413  fe2.getOpPtrVector().push_back(
414  new MyOp2(my_split, UMDataOp::OPROWCOL, ContactDataOp::FACESLAVESLAVE));
415  CHKERR m_field.loop_finite_elements("TEST_PROBLEM", "TEST_FE2", fe2);
416  }
417  CATCH_ERRORS;
418 
420 
421  return 0;
422 }
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
CoreTmp< 0 > Core
Definition: Core.hpp:1129
scalar or vector of scalars describe (no true field)
Definition: definitions.h:176
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:485
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:60
MoFEMErrorCode splitSides(const EntityHandle meshset, const BitRefLevel &bit, const int msId, const CubitBCType cubit_bc_type, const bool add_interface_entities, const bool recursive=false, int verb=QUIET)
Split nodes and other entities of tetrahedra on both sides of the interface and insert flat prisms in...
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.
char mesh_file_name[255]
Projection of edge entities with one mid-node on hierarchical basis.
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
MoFEMErrorCode getInterface(const MOFEMuuid &uuid, IFACE *&iface) const
Get interface by uuid and return reference to pointer of interface.
ForcesAndSourcesCore::UserDataOperator UserDataOperator
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:152
Create interface from given surface and insert flat prisms in-between.
tee_device< std::ostream, std::ofstream > TeeDevice
virtual MoFEMErrorCode loop_dofs(const Problem *problem_ptr, const std::string &field_name, RowColData rc, DofMethod &method, int lower_rank, int upper_rank, int verb=DEFAULT_VERBOSITY)=0
Make a loop over dofs.
virtual MoFEMErrorCode modify_problem_ref_level_add_bit(const std::string &name_problem, const BitRefLevel &bit)=0
add ref level to problem
MoFEMErrorCode getSides(const int msId, const CubitBCType cubit_bc_type, const BitRefLevel mesh_bit_level, const bool recursive, int verb=QUIET)
Store tetrahedra from each side of the interface separately in two child meshsets of the parent meshs...
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
#define _IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet having a particular BC meshset in a moFEM field.
#define CHKERR
Inline error check.
Definition: definitions.h:604
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
virtual EntityHandle get_field_meshset(const std::string name) const =0
get field meshset
constexpr int order
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:292
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
ContactPrism finite elementUser is implementing own operator at Gauss points level,...
virtual MoFEMErrorCode add_problem(const std::string &name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add problem.
#define _IT_CUBITMESHSETS_FOR_LOOP_(MESHSET_MANAGER, IT)
Iterator that loops over all the Cubit MeshSets in a moFEM field.
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1943
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....
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
Core (interface) class.
Definition: Core.hpp:77
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:415
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:441
continuous field
Definition: definitions.h:177
MoFEMErrorCode partitionSimpleProblem(const std::string name, int verb=VERBOSE)
partition problem dofs
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:100
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...
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)

Variable Documentation

◆ help

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