v0.9.0
forces_and_sources_testing_edge_element.cpp
Go to the documentation of this file.
1 /** \file forces_and_sources_testing_edge_element.cpp
2  * \example forces_and_sources_testing_edge_element.cpp
3  * \brief Testing edge elements
4  * nodesets.
5  *
6 */
7 
8 /* This file is part of MoFEM.
9  * MoFEM is free software: you can redistribute it and/or modify it under
10  * the terms of the GNU Lesser General Public License as published by the
11  * Free Software Foundation, either version 3 of the License, or (at your
12  * option) any later version.
13  *
14  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
15  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
16  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
17  * License for more details.
18  *
19  * You should have received a copy of the GNU Lesser General Public
20  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
21 
22 #include <MoFEM.hpp>
23 
24 namespace bio = boost::iostreams;
25 using bio::stream;
26 using bio::tee_device;
27 
28 using namespace MoFEM;
29 
30 static char help[] = "...\n\n";
31 
32 int main(int argc, char *argv[]) {
33 
34  MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
35 
36  try {
37 
38  moab::Core mb_instance;
39  moab::Interface &moab = mb_instance;
40  int rank;
41  MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
42 
43  PetscBool flg = PETSC_TRUE;
44  char mesh_file_name[255];
45 #if PETSC_VERSION_GE(3, 6, 4)
46  CHKERR PetscOptionsGetString(PETSC_NULL, "", "-my_file", mesh_file_name,
47  255, &flg);
48 
49 #else
50  CHKERR PetscOptionsGetString(PETSC_NULL, PETSC_NULL, "-my_file",
51  mesh_file_name, 255, &flg);
52 
53 #endif
54 
55  if (flg != PETSC_TRUE) {
56  SETERRQ(PETSC_COMM_SELF, 1, "*** ERROR -my_file (MESH FILE NEEDED)");
57  }
58 
59  // Create MoFEM database
60  MoFEM::Core core(moab);
61  MoFEM::Interface &m_field = core;
62 
63  ParallelComm *pcomm = ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
64  if (pcomm == NULL)
65  pcomm = new ParallelComm(&moab, PETSC_COMM_WORLD);
66 
67  const char *option;
68  option = ""; //"PARALLEL=BCAST;";//;DEBUG_IO";
69  CHKERR moab.load_file(mesh_file_name, 0, option);
70 
71  // set ebturues bit level
72  BitRefLevel bit_level0;
73  bit_level0.set(0);
74  CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevelByDim(
75  0, 3, bit_level0);
76 
77  enum bases { AINSWORTH, BERNSTEIN_BEZIER, LASBASETOP };
78  const char *list_bases[] = {"ainsworth", "bernstein_bezier"};
79  PetscInt choice_base_value = AINSWORTH;
80  CHKERR PetscOptionsGetEList(PETSC_NULL, NULL, "-base", list_bases,
81  LASBASETOP, &choice_base_value, &flg);
83  if (choice_base_value == AINSWORTH)
85  else if (choice_base_value == BERNSTEIN_BEZIER)
87 
88  // Fields
89  CHKERR m_field.add_field("FIELD1", H1, base, 3);
90  CHKERR m_field.add_field("MESH_NODE_POSITIONS", H1, AINSWORTH_LEGENDRE_BASE,
91  3);
92 
93  // FE
94  CHKERR m_field.add_finite_element("TEST_FE");
95 
96 
97  // Define rows/cols and element data
98  auto add_field_to_fe = [&m_field](const std::string field_name) {
100  CHKERR m_field.modify_finite_element_add_field_row("TEST_FE", field_name);
101  CHKERR m_field.modify_finite_element_add_field_col("TEST_FE", field_name);
103  field_name);
105  };
106  CHKERR add_field_to_fe("FIELD1");
108  "MESH_NODE_POSITIONS");
109 
110  // Problem
111  CHKERR m_field.add_problem("TEST_PROBLEM");
112 
113  // set finite elements for problem
114  CHKERR m_field.modify_problem_add_finite_element("TEST_PROBLEM", "TEST_FE");
115  // set refinement level for problem
116  CHKERR m_field.modify_problem_ref_level_add_bit("TEST_PROBLEM", bit_level0);
117 
118 
119  // meshset consisting all entities in mesh
120  EntityHandle root_set = moab.get_root_set();
121  // add entities to field
122  CHKERR m_field.add_ents_to_field_by_type(root_set, MBEDGE, "FIELD1");
123  CHKERR m_field.add_ents_to_field_by_type(root_set, MBEDGE,
124  "MESH_NODE_POSITIONS");
125  CHKERR m_field.set_field_order(root_set, MBVERTEX, "MESH_NODE_POSITIONS",
126  1);
127  CHKERR m_field.set_field_order(root_set, MBEDGE, "MESH_NODE_POSITIONS", 1);
128  CHKERR m_field.set_field_order(root_set, MBTRI, "MESH_NODE_POSITIONS", 1);
129  CHKERR m_field.set_field_order(root_set, MBTET, "MESH_NODE_POSITIONS", 1);
130 
131 
132  // add entities to finite element
133  Range tets;
134  CHKERR moab.get_entities_by_type(0, MBTET, tets, false);
135  Skinner skin(&m_field.get_moab());
136  Range tets_skin;
137  CHKERR skin.find_skin(0, tets, false, tets_skin);
138  Range tets_skin_edges;
139  CHKERR moab.get_adjacencies(tets_skin, 1, false, tets_skin_edges,
140  moab::Interface::UNION);
141  CHKERR m_field.add_ents_to_finite_element_by_type(tets_skin_edges, MBEDGE,
142  "TEST_FE");
143 
144  // set app. order
145  // see Hierarchic Finite Element Bases on Unstructured Tetrahedral Meshes
146  // (Mark Ainsworth & Joe Coyle)
147  int order = 3;
148  CHKERR m_field.set_field_order(root_set, MBEDGE, "FIELD1", order);
149  CHKERR m_field.set_field_order(root_set, MBVERTEX, "FIELD1", 1);
150 
151  /****/
152  // build database
153  // build field
154  CHKERR m_field.build_fields();
155 
156  // set FIELD1 from positions of 10 node tets
157  if (base == AINSWORTH_LEGENDRE_BASE) {
158  Projection10NodeCoordsOnField ent_method(m_field, "FIELD1");
159  CHKERR m_field.loop_dofs("FIELD1", ent_method);
160 
161  }
162  {
163  Projection10NodeCoordsOnField ent_method(m_field, "MESH_NODE_POSITIONS");
164  CHKERR m_field.loop_dofs("MESH_NODE_POSITIONS", ent_method);
165 
166  }
167  // build finite elemnts
168  CHKERR m_field.build_finite_elements();
169 
170  // build adjacencies
171  CHKERR m_field.build_adjacencies(bit_level0);
172 
173  // build problem
174  ProblemsManager *prb_mng_ptr;
175  CHKERR m_field.getInterface(prb_mng_ptr);
176 
177  CHKERR prb_mng_ptr->buildProblem("TEST_PROBLEM", true);
178 
179  // mesh partitioning
180  // partition
181  CHKERR prb_mng_ptr->partitionSimpleProblem("TEST_PROBLEM");
182  CHKERR prb_mng_ptr->partitionFiniteElements("TEST_PROBLEM");
183 
184  // what are ghost nodes, see Petsc Manual
185  CHKERR prb_mng_ptr->partitionGhostDofs("TEST_PROBLEM");
186 
187  EdgeElementForcesAndSourcesCore fe1(m_field);
188 
189  typedef tee_device<std::ostream, std::ofstream> TeeDevice;
190  typedef stream<TeeDevice> TeeStream;
191 
192  std::ofstream ofs(
193  (std ::string("forces_and_sources_testing_edge_element_") +
194  ApproximationBaseNames[base] + ".txt")
195  .c_str());
196  TeeDevice my_tee(std::cout, ofs);
197  TeeStream my_split(my_tee);
198 
200 
201  TeeStream &my_split;
202  MyOp(TeeStream &_my_split, const char type)
204  "FIELD1", type),
205  my_split(_my_split) {}
206 
207  MoFEMErrorCode doWork(int side, EntityType type,
210 
211  my_split << "NH1" << std::endl;
212  my_split << "side: " << side << " type: " << type << std::endl;
213  my_split << "data: " << data << std::endl;
214  my_split << "coords: " << std::setprecision(3) << getCoords()
215  << std::endl;
216  my_split << "getCoordsAtGaussPts: " << std::setprecision(3)
217  << getCoordsAtGaussPts() << std::endl;
218  my_split << "length: " << std::setprecision(3) << getLength()
219  << std::endl;
220  my_split << "direction: " << std::setprecision(3) << getDirection()
221  << std::endl;
222 
223  int nb_gauss_pts = data.getN().size1();
224  for (int gg = 0; gg < nb_gauss_pts; gg++) {
225  my_split << "tangent " << gg << " " << getTangetAtGaussPts()
226  << std::endl;
227  }
228 
230  }
231 
232  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
233  EntityType col_type,
237  my_split << "ROW NH1NH1" << std::endl;
238  my_split << "row side: " << row_side << " row_type: " << row_type
239  << std::endl;
240  my_split << row_data << std::endl;
241  my_split << "COL NH1NH1" << std::endl;
242  my_split << "col side: " << col_side << " col_type: " << col_type
243  << std::endl;
244  my_split << col_data << std::endl;
246  }
247  };
248 
249  fe1.getOpPtrVector().push_back(
251  fe1.getOpPtrVector().push_back(
253 
254  CHKERR m_field.loop_finite_elements("TEST_PROBLEM", "TEST_FE", fe1);
255  }
256  CATCH_ERRORS;
257 
259 
260  return 0;
261 }
int main(int argc, char *argv[])
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
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
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:477
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...
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
char mesh_file_name[255]
Projection of edge entities with one mid-node on hierarchical basis.
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
static const char *const ApproximationBaseNames[]
Definition: definitions.h:158
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.
FieldApproximationBase
approximation base
Definition: definitions.h:144
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
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.
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:66
Managing BitRefLevels.
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:146
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
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
boost::ptr_vector< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.
#define CHKERR
Inline error check.
Definition: definitions.h:596
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:285
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
virtual MoFEMErrorCode add_problem(const std::string &name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add problem.
Data on single entity (This is passed as argument to DataOperator::doWork)
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
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:433
continuous field
Definition: definitions.h:171
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
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)