v0.14.0
Loading...
Searching...
No Matches
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
9
10#include <MoFEM.hpp>
11
12namespace bio = boost::iostreams;
13using bio::stream;
14using bio::tee_device;
15
16using namespace MoFEM;
17
18static char help[] = "...\n\n";
19
20int main(int argc, char *argv[]) {
21
22 MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
23
24 try {
25
26 moab::Core mb_instance;
27 moab::Interface &moab = mb_instance;
28 int rank;
29 MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
30
31 PetscBool flg = PETSC_TRUE;
32 char mesh_file_name[255];
33#if PETSC_VERSION_GE(3, 6, 4)
34 CHKERR PetscOptionsGetString(PETSC_NULL, "", "-my_file", mesh_file_name,
35 255, &flg);
36
37#else
38 CHKERR PetscOptionsGetString(PETSC_NULL, PETSC_NULL, "-my_file",
39 mesh_file_name, 255, &flg);
40
41#endif
42
43 if (flg != PETSC_TRUE) {
44 SETERRQ(PETSC_COMM_SELF, 1, "*** ERROR -my_file (MESH FILE NEEDED)");
45 }
46
47 // Create MoFEM database
48 MoFEM::Core core(moab);
49 MoFEM::Interface &m_field = core;
50
51 const char *option;
52 option = ""; ;
53 CHKERR moab.load_file(mesh_file_name, 0, option);
54
55 // set ebturues bit level
56 BitRefLevel bit_level0;
57 bit_level0.set(0);
58 CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevelByDim(
59 0, 3, bit_level0);
60
61 enum bases { AINSWORTH, BERNSTEIN_BEZIER, LASBASETOP };
62 const char *list_bases[] = {"ainsworth", "bernstein_bezier"};
63 PetscInt choice_base_value = AINSWORTH;
64 CHKERR PetscOptionsGetEList(PETSC_NULL, NULL, "-base", list_bases,
65 LASBASETOP, &choice_base_value, &flg);
67 if (choice_base_value == AINSWORTH)
69 else if (choice_base_value == BERNSTEIN_BEZIER)
71
72 // Fields
73 CHKERR m_field.add_field("FIELD1", H1, base, 3);
74 CHKERR m_field.add_field("MESH_NODE_POSITIONS", H1, AINSWORTH_LEGENDRE_BASE,
75 3);
76
77 // FE
78 CHKERR m_field.add_finite_element("TEST_FE");
79
80
81 // Define rows/cols and element data
82 auto add_field_to_fe = [&m_field](const std::string field_name) {
89 };
90 CHKERR add_field_to_fe("FIELD1");
92 "MESH_NODE_POSITIONS");
93
94 // Problem
95 CHKERR m_field.add_problem("TEST_PROBLEM");
96
97 // set finite elements for problem
98 CHKERR m_field.modify_problem_add_finite_element("TEST_PROBLEM", "TEST_FE");
99 // set refinement level for problem
100 CHKERR m_field.modify_problem_ref_level_add_bit("TEST_PROBLEM", bit_level0);
101
102
103 // meshset consisting all entities in mesh
104 EntityHandle root_set = moab.get_root_set();
105 // add entities to field
106 CHKERR m_field.add_ents_to_field_by_type(root_set, MBEDGE, "FIELD1");
107 CHKERR m_field.add_ents_to_field_by_type(root_set, MBEDGE,
108 "MESH_NODE_POSITIONS");
109 CHKERR m_field.set_field_order(root_set, MBVERTEX, "MESH_NODE_POSITIONS",
110 1);
111 CHKERR m_field.set_field_order(root_set, MBEDGE, "MESH_NODE_POSITIONS", 1);
112 CHKERR m_field.set_field_order(root_set, MBTRI, "MESH_NODE_POSITIONS", 1);
113 CHKERR m_field.set_field_order(root_set, MBTET, "MESH_NODE_POSITIONS", 1);
114
115
116 // add entities to finite element
117 Range tets;
118 CHKERR moab.get_entities_by_type(0, MBTET, tets, false);
119 Skinner skin(&m_field.get_moab());
120 Range tets_skin;
121 CHKERR skin.find_skin(0, tets, false, tets_skin);
122 Range tets_skin_edges;
123 CHKERR moab.get_adjacencies(tets_skin, 1, false, tets_skin_edges,
124 moab::Interface::UNION);
125 CHKERR m_field.add_ents_to_finite_element_by_type(tets_skin_edges, MBEDGE,
126 "TEST_FE");
127
128 // set app. order
129 // see Hierarchic Finite Element Bases on Unstructured Tetrahedral Meshes
130 // (Mark Ainsworth & Joe Coyle)
131 int order = 3;
132 CHKERR m_field.set_field_order(root_set, MBEDGE, "FIELD1", order);
133 CHKERR m_field.set_field_order(root_set, MBVERTEX, "FIELD1", 1);
134
135 /****/
136 // build database
137 // build field
138 CHKERR m_field.build_fields();
139
140 // set FIELD1 from positions of 10 node tets
141 if (base == AINSWORTH_LEGENDRE_BASE) {
142 Projection10NodeCoordsOnField ent_method(m_field, "FIELD1");
143 CHKERR m_field.loop_dofs("FIELD1", ent_method);
144
145 }
146 {
147 Projection10NodeCoordsOnField ent_method(m_field, "MESH_NODE_POSITIONS");
148 CHKERR m_field.loop_dofs("MESH_NODE_POSITIONS", ent_method);
149
150 }
151 // build finite elemnts
153
154 // build adjacencies
155 CHKERR m_field.build_adjacencies(bit_level0);
156
157 // build problem
158 ProblemsManager *prb_mng_ptr;
159 CHKERR m_field.getInterface(prb_mng_ptr);
160
161 CHKERR prb_mng_ptr->buildProblem("TEST_PROBLEM", true);
162
163 // mesh partitioning
164 // partition
165 CHKERR prb_mng_ptr->partitionSimpleProblem("TEST_PROBLEM");
166 CHKERR prb_mng_ptr->partitionFiniteElements("TEST_PROBLEM");
167
168 // what are ghost nodes, see Petsc Manual
169 CHKERR prb_mng_ptr->partitionGhostDofs("TEST_PROBLEM");
170
172
173 typedef tee_device<std::ostream, std::ofstream> TeeDevice;
174 typedef stream<TeeDevice> TeeStream;
175
176 std::ofstream ofs(
177 (std ::string("forces_and_sources_testing_edge_element_") +
178 ApproximationBaseNames[base] + ".txt")
179 .c_str());
180 TeeDevice my_tee(std::cout, ofs);
181 TeeStream my_split(my_tee);
182
183 struct MyOp : public EdgeElementForcesAndSourcesCore::UserDataOperator {
184
185 TeeStream &my_split;
186 MyOp(TeeStream &_my_split, const char type)
187 : EdgeElementForcesAndSourcesCore::UserDataOperator("FIELD1",
188 "FIELD1", type),
189 my_split(_my_split) {}
190
191 MoFEMErrorCode doWork(int side, EntityType type,
194
195 my_split << "NH1" << std::endl;
196 my_split << "side: " << side << " type: " << type << std::endl;
197 my_split << "data: " << data << std::endl;
198 my_split << "coords: " << std::setprecision(3) << getCoords()
199 << std::endl;
200 my_split << "getCoordsAtGaussPts: " << std::setprecision(3)
201 << getCoordsAtGaussPts() << std::endl;
202 my_split << "length: " << std::setprecision(3) << getLength()
203 << std::endl;
204 my_split << "direction: " << std::setprecision(3) << getDirection()
205 << std::endl;
206
207 int nb_gauss_pts = data.getN().size1();
208 for (int gg = 0; gg < nb_gauss_pts; gg++) {
209 my_split << "tangent " << gg << " " << getTangentAtGaussPts()
210 << std::endl;
211 }
212
214 }
215
216 MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
217 EntityType col_type,
219 EntitiesFieldData::EntData &col_data) {
221 my_split << "ROW NH1NH1" << std::endl;
222 my_split << "row side: " << row_side << " row_type: " << row_type
223 << std::endl;
224 my_split << row_data << std::endl;
225 my_split << "COL NH1NH1" << std::endl;
226 my_split << "col side: " << col_side << " col_type: " << col_type
227 << std::endl;
228 my_split << col_data << std::endl;
230 }
231 };
232
233 fe1.getOpPtrVector().push_back(
234 new OpGetHOTangentsOnEdge("MESH_NODE_POSITIONS"));
235 fe1.getOpPtrVector().push_back(
236 new MyOp(my_split, ForcesAndSourcesCore::UserDataOperator::OPROW));
237 fe1.getOpPtrVector().push_back(
238 new MyOp(my_split, ForcesAndSourcesCore::UserDataOperator::OPROWCOL));
239
240 CHKERR m_field.loop_finite_elements("TEST_PROBLEM", "TEST_FE", fe1);
241 }
243
245
246 return 0;
247}
int main()
#define CATCH_ERRORS
Catch errors.
FieldApproximationBase
approximation base
Definition definitions.h:58
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition definitions.h:60
@ NOBASE
Definition definitions.h:59
@ AINSWORTH_BERNSTEIN_BEZIER_BASE
Definition definitions.h:64
@ H1
continuous field
Definition definitions.h:85
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
static const char *const ApproximationBaseNames[]
Definition definitions.h:72
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
constexpr int order
tee_device< std::ostream, std::ofstream > TeeDevice
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
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
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string name_filed)=0
set finite element field data
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
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
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.
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.
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 loop_finite_elements(const std::string problem_name, const std::string &fe_name, FEMethod &method, boost::shared_ptr< NumeredEntFiniteElement_multiIndex > fe_ptr=nullptr, MoFEMTypes bh=MF_EXIST, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr(), int verb=DEFAULT_VERBOSITY)=0
Make a loop over finite elements.
MoFEMErrorCode partitionGhostDofs(const std::string name, int verb=VERBOSE)
determine ghost nodes
MoFEMErrorCode partitionSimpleProblem(const std::string name, int verb=VERBOSE)
partition problem dofs
MoFEMErrorCode buildProblem(const std::string name, const bool square_matrix, int verb=VERBOSE)
build problem data structures
MoFEMErrorCode partitionFiniteElements(const std::string name, bool part_from_moab=false, int low_proc=-1, int hi_proc=-1, int verb=VERBOSE)
partition finite elements
virtual MoFEMErrorCode add_problem(const std::string &name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add problem.
virtual MoFEMErrorCode modify_problem_ref_level_add_bit(const std::string &name_problem, const BitRefLevel &bit)=0
add ref level to problem
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
char mesh_file_name[255]
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition Types.hpp:40
implementation of Data Operators for Forces and Sources
Definition Common.hpp:10
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
constexpr auto field_name
Managing BitRefLevels.
virtual moab::Interface & get_moab()=0
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
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.
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.
Data on single entity (This is passed as argument to DataOperator::doWork)
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
MatrixDouble & getCoordsAtGaussPts()
Gauss points and weight, matrix (nb. of points x 3)
boost::ptr_deque< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.
Calculate tangent vector on edge form HO geometry approximation.
Problem manager is used to build and partition problems.
Projection of edge entities with one mid-node on hierarchical basis.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Operator used to check consistency between local coordinates and global cooridnates for integrated po...
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
MyOp(std::array< double, 12 > &eval_points)