v0.13.1
Loading...
Searching...
No Matches
h_adaptive_transport.cpp
Go to the documentation of this file.
1/** \file h_adaptive_transport.cpp
2\brief Example implementation of transport problem using mixed formulation
3
4\todo Should be implemented and tested problem from this article
5Demkowicz, Leszek, and Jayadeep Gopalakrishnan. "Analysis of the DPG method for
6the Poisson equation." SIAM Journal on Numerical Analysis 49.5 (2011):
71788-1809.
8
9\ingroup mofem_mix_transport_elem
10*/
11
12
13
16
17using namespace MoFEM;
18using namespace MixTransport;
19
20static char help[] = "-my_file input file"
21 "-my_order order of approximation"
22 "-nb_levels number of refinements levels"
23 "\n\n";
24
25/**
26 * Data structure to pass information between function evaluating boundary
27 * values and fluxes and generic data structures for boundary conditions on
28 * meshsets.
29 */
30struct BcFluxData {
32 double fLux;
33};
34typedef map<int, BcFluxData> BcFluxMap;
35
36/** \brief Application of mix transport data structure
37 *
38 * MixTransportElement is a class collecting functions, operators and
39 * data for mix implementation of transport element. See there to
40 * learn how elements are created or how operators look like.
41 *
42 * Some methods in MixTransportElement are abstract, f.e. user need to
43 * implement own source therm.
44
45 * \ingroup mofem_mix_transport_elem
46 */
47struct MyTransport : public MixTransportElement {
48
51 double lastFlux;
52
53 MyTransport(MoFEM::Interface &m_field, BcFluxMap &bc_flux_map)
54 : MixTransportElement(m_field), bcFluxMap(bc_flux_map), lastEnt(0),
55 lastFlux(0) {}
56
57 /**
58 * \brief set source term
59 * @param ent handle to entity on which function is evaluated
60 * @param x coord
61 * @param y coord
62 * @param z coord
63 * @param flux reference to source term set by function
64 * @return error code
65 */
66 MoFEMErrorCode getSource(EntityHandle ent, const double x, const double y,
67 const double z, double &flux) {
69 flux = 0;
71 }
72
73 /**
74 * \brief natural (Dirihlet) boundary conditions (set values)
75 * @param ent handle to finite element entity
76 * @param x coord
77 * @param y coord
78 * @param z coord
79 * @param value reference to value set by function
80 * @return error code
81 */
82 MoFEMErrorCode getBcOnValues(const EntityHandle ent, const double x,
83 const double y, const double z, double &value) {
85 value = 0;
87 }
88
89 /**
90 * \brief essential (Neumann) boundary condition (set fluxes)
91 * @param ent handle to finite element entity
92 * @param x coord
93 * @param y coord
94 * @param z coord
95 * @param flux reference to flux which is set by function
96 * @return [description]
97 */
98 MoFEMErrorCode getBcOnFluxes(const EntityHandle ent, const double x,
99 const double y, const double z, double &flux) {
101 if (lastEnt == ent) {
102 flux = lastFlux;
103 } else {
104 flux = 0;
105 for (BcFluxMap::iterator mit = bcFluxMap.begin(); mit != bcFluxMap.end();
106 mit++) {
107 Range &tris = mit->second.eNts;
108 if (tris.find(ent) != tris.end()) {
109 flux = mit->second.fLux;
110 }
111 }
112 lastEnt = ent;
113 lastFlux = flux;
114 }
116 }
117
118 /**
119 * \brief set-up boundary conditions
120 * @param ref_level mesh refinement level
121
122 \note It is assumed that user would like to something non-standard with
123 boundary conditions, have a own type of data structures to pass to functions
124 calculating values and fluxes on boundary. For example BcFluxMap. That way
125 this function is implemented here not in generic class MixTransportElement.
126
127 * @return error code
128 */
131 Range tets;
132 CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
133 ref_level, BitRefLevel().set(), MBTET, tets);
134 Skinner skin(&mField.get_moab());
135 Range skin_faces; // skin faces from 3d ents
136 CHKERR skin.find_skin(0, tets, false, skin_faces);
137 // note: what is essential (dirichlet) is natural (neumann) for mix-FE
138 // compared to classical FE
139 Range natural_bc;
141 mField, NODESET | TEMPERATURESET, it)) {
142 Range tris;
143 CHKERR it->getMeshsetIdEntitiesByDimension(mField.get_moab(), 2, tris,
144 true);
145 natural_bc.insert(tris.begin(), tris.end());
146 }
148 mField, SIDESET | HEATFLUXSET, it)) {
149 HeatFluxCubitBcData mydata;
150 CHKERR it->getBcDataStructure(mydata);
151 if (mydata.data.flag1 == 1) {
152 Range tris;
153 CHKERR it->getMeshsetIdEntitiesByDimension(mField.get_moab(), 2, tris,
154 true);
155 bcFluxMap[it->getMeshsetId()].eNts = tris;
156 bcFluxMap[it->getMeshsetId()].fLux = mydata.data.value1;
157 // cerr << bcFluxMap[it->getMeshsetId()].eNts << endl;
158 // cerr << bcFluxMap[it->getMeshsetId()].fLux << endl;
159 }
160 }
161 Range essential_bc = subtract(skin_faces, natural_bc);
162 Range bit_tris;
163 CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
164 ref_level, BitRefLevel().set(), MBTRI, bit_tris);
165 essential_bc = intersect(bit_tris, essential_bc);
166 natural_bc = intersect(bit_tris, natural_bc);
168 "MIX_BCFLUX");
170 "MIX_BCVALUE");
171 // CHKERR
172 // mField.add_ents_to_finite_element_by_type(skin_faces,MBTRI,"MIX_BCVALUE");
174 }
175
176 /**
177 * \brief Refine mesh
178 * @param ufe general data structure
179 * @param nb_levels number of refinement levels
180 * @param order set order of approximation
181 * @return errpr code
182
183 Refinement of could result in distorted mesh, for example, imagine when you
184 have two levels of non-uniform refinement. Some tetrahedra on the mesh at
185 first refinement instance are only refined by splitting subset of edges on
186 it. Then refined child tetrahedra usually will have worse quality than
187 quality of parent element. Refining such element in subsequent mesh
188 refinement, potentially will deteriorate elements quality even worse. To
189 prevent that adding new refinement level, recreate whole hierarchy of meshes.
190
191 Note on subsequent improvement could include refinement of
192 tetrahedra from different levels, including initial mesh. So refinement two
193 could split elements created during refinement one and also split elements
194 from an initial mesh.
195
196 That adding the new refinement level creates refinement hierarchy of meshes
197 from a scratch, not adding to existing one.
198
199 Entities from previous hierarchy are used in that process, but bit levels on
200 those entities are squashed.
201
202 */
204 const int order) {
205 MeshRefinement *refine_ptr;
207 // get refined edges having child vertex
208 auto ref_ents_ptr = mField.get_ref_ents();
209 typedef RefEntity_multiIndex::index<
210 Composite_EntType_and_ParentEntType_mi_tag>::type RefEntsByComposite;
211 const RefEntsByComposite &ref_ents =
213 RefEntsByComposite::iterator rit, hi_rit;
214 rit = ref_ents.lower_bound(boost::make_tuple(MBVERTEX, MBEDGE));
215 hi_rit = ref_ents.upper_bound(boost::make_tuple(MBVERTEX, MBEDGE));
216 Range refined_edges;
217 // thist loop is over vertices which parent is edge
218 for (; rit != hi_rit; rit++) {
219 refined_edges.insert((*rit)->getParentEnt()); // get parent edge
220 }
221 // get tets which has large error
222 Range tets_to_refine;
223 const double max_error = ufe.errorMap.rbegin()->first;
224 // int size = ((double)5/6)*ufe.errorMap.size();
225 for (map<double, EntityHandle>::iterator mit = ufe.errorMap.begin();
226 mit != ufe.errorMap.end(); mit++) {
227 // cerr << mit->first << " " << mit->second << endl;
228 // if((size--)>0) continue;
229 if (mit->first < 0.25 * max_error)
230 continue;
231 tets_to_refine.insert(mit->second);
232 }
233 Range tets_to_refine_edges;
234 CHKERR mField.get_moab().get_adjacencies(
235 tets_to_refine, 1, false, tets_to_refine_edges, moab::Interface::UNION);
236 refined_edges.merge(tets_to_refine_edges);
237 CHKERR mField.getInterface(refine_ptr);
238 for (int ll = 0; ll != nb_levels; ll++) {
239 Range edges;
240 CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
241 BitRefLevel().set(ll), BitRefLevel().set(), MBEDGE, edges);
242 edges = intersect(edges, refined_edges);
243 // add edges to refine at current level edges (some of the where refined
244 // before)
246 edges, BitRefLevel().set(ll + 1));
247 // get tets at current level
248 Range tets;
249 CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
250 BitRefLevel().set(ll), BitRefLevel().set(), MBTET, tets);
251 CHKERR refine_ptr->refineTets(tets, BitRefLevel().set(ll + 1));
253 }
254
255 // update fields and elements
256 EntityHandle ref_meshset;
257 CHKERR mField.get_moab().create_meshset(MESHSET_SET, ref_meshset);
258 {
259 // cerr << BitRefLevel().set(nb_levels) << endl;
260 CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
261 BitRefLevel().set(nb_levels), BitRefLevel().set(), MBTET,
262 ref_meshset);
263
264 Range ref_tets;
265 CHKERR mField.get_moab().get_entities_by_type(ref_meshset, MBTET,
266 ref_tets);
267
268 // add entities to field
269 CHKERR mField.add_ents_to_field_by_type(ref_meshset, MBTET, "FLUXES");
270 CHKERR mField.add_ents_to_field_by_type(ref_meshset, MBTET, "VALUES");
271 CHKERR mField.set_field_order(0, MBTET, "FLUXES", order + 1);
272 CHKERR mField.set_field_order(0, MBTRI, "FLUXES", order + 1);
273 CHKERR mField.set_field_order(0, MBTET, "VALUES", order);
274
275 // add entities to skeleton
276 Range ref_tris;
277 CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
278 BitRefLevel().set(nb_levels), BitRefLevel().set(), MBTRI, ref_tris);
280 "MIX_SKELETON");
281
282 // add entities to finite elements
285 Mat_Thermal temp_data;
286 CHKERR it->getAttributeDataStructure(temp_data);
287 setOfBlocks[it->getMeshsetId()].cOnductivity =
288 temp_data.data.Conductivity;
289 setOfBlocks[it->getMeshsetId()].cApacity = temp_data.data.HeatCapacity;
290 CHKERR mField.get_moab().get_entities_by_type(
291 it->meshset, MBTET, setOfBlocks[it->getMeshsetId()].tEts, true);
292 setOfBlocks[it->getMeshsetId()].tEts =
293 intersect(ref_tets, setOfBlocks[it->getMeshsetId()].tEts);
295 setOfBlocks[it->getMeshsetId()].tEts, MBTET, "MIX");
296 }
297 }
298 CHKERR mField.get_moab().delete_entities(&ref_meshset, 1);
300 }
301
302 /**
303 * \brief Squash bits of entities
304
305 Information about hierarchy of meshsets is lost, but entities are not deleted
306 from the mesh. After squash entities bits, new hierarchy can be created.
307
308 * @return error code
309 */
312 BitRefLevel all_but_0;
313 all_but_0.set(0);
314 all_but_0.flip();
315 BitRefLevel garbage_bit;
316 garbage_bit.set(BITREFLEVEL_SIZE - 1); // Garbage level
317 auto ref_ents_ptr = mField.get_ref_ents();
318 RefEntity_multiIndex::iterator mit = ref_ents_ptr->begin();
319 for (; mit != ref_ents_ptr->end(); mit++) {
320 if (mit->get()->getEntType() == MBENTITYSET)
321 continue;
322 BitRefLevel bit = mit->get()->getBitRefLevel();
323 if ((all_but_0 & bit) == bit) {
324 *(const_cast<RefEntity *>(mit->get())->getBitRefLevelPtr()) =
325 garbage_bit;
326 } else {
327 *(const_cast<RefEntity *>(mit->get())->getBitRefLevelPtr()) =
328 BitRefLevel().set(0);
329 }
330 }
332 }
333
334 /**
335 * \brief update meshsets with new entities after mesh refinement
336 * @param nb_levels nb_levels
337 * @param order appropriate order
338 * @return error code
339 */
341 BitRefLevel ref_level;
343 ref_level.set(nb_levels);
345 EntityHandle meshset = it->meshset;
347 ->updateMeshsetByEntitiesChildren(meshset, ref_level, meshset,
348 MBMAXTYPE, true);
349 }
351 }
352};
353
354int main(int argc, char *argv[]) {
355
356 const string default_options = "-ksp_type fgmres \n"
357 "-pc_type lu \n"
358 "-pc_factor_mat_solver_type mumps \n"
359 "-mat_mumps_icntl_20 0 \n"
360 "-ksp_monitor\n";
361
362 string param_file = "param_file.petsc";
363 if (!static_cast<bool>(ifstream(param_file))) {
364 std::ofstream file(param_file.c_str(), std::ios::ate);
365 if (file.is_open()) {
366 file << default_options;
367 file.close();
368 }
369 }
370
371 MoFEM::Core::Initialize(&argc, &argv, param_file.c_str(), help);
372
373 try {
374
375 moab::Core mb_instance;
376 moab::Interface &moab = mb_instance;
377 int rank;
378 MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
379
380 // get file name form command line
381 PetscBool flg = PETSC_TRUE;
382 char mesh_file_name[255];
383 CHKERR PetscOptionsGetString(PETSC_NULL, PETSC_NULL, "-my_file",
384 mesh_file_name, 255, &flg);
385 if (flg != PETSC_TRUE) {
386 SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
387 "*** ERROR -my_file (MESH FILE NEEDED)");
388 }
389
390 const char *option;
391 option = "";
392 CHKERR moab.load_file(mesh_file_name, 0, option);
393
394 // Create mofem interface
395 MoFEM::Core core(moab);
396 MoFEM::Interface &m_field = core;
397
398 // Add meshsets with material and boundary conditions
399 MeshsetsManager *meshsets_manager_ptr;
400 CHKERR m_field.getInterface(meshsets_manager_ptr);
401 CHKERR meshsets_manager_ptr->setMeshsetFromFile();
402
403 PetscPrintf(PETSC_COMM_WORLD,
404 "Read meshsets add added meshsets for bc.cfg\n");
405 for (_IT_CUBITMESHSETS_FOR_LOOP_(m_field, it)) {
406 PetscPrintf(PETSC_COMM_WORLD, "%s",
407 static_cast<std::ostringstream &>(
408 std::ostringstream().seekp(0) << *it << endl)
409 .str()
410 .c_str());
411 cerr << *it << endl;
412 }
413
414 // set entities bit level
415 BitRefLevel ref_level;
416 ref_level.set(0);
417 CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevelByDim(
418 0, 3, ref_level);
419
420 // set app. order
421 // see Hierarchic Finite Element Bases on Unstructured Tetrahedral Meshes
422 // (Mark Ainsworth & Joe Coyle)
423 PetscInt order;
424 CHKERR PetscOptionsGetInt(PETSC_NULL, PETSC_NULL, "-my_order", &order,
425 &flg);
426 if (flg != PETSC_TRUE) {
427 order = 0;
428 }
429
430 // finite elements
431
432 BcFluxMap bc_flux_map;
433 MyTransport ufe(m_field, bc_flux_map);
434
435 // Initially calculate problem on coarse mesh
436
437 CHKERR ufe.addFields("VALUES", "FLUXES", order);
438 CHKERR ufe.addFiniteElements("FLUXES", "VALUES");
439 // Set boundary conditions
440 CHKERR ufe.addBoundaryElements(ref_level);
441 CHKERR ufe.buildProblem(ref_level);
445 CHKERR ufe.evaluateError();
447 CHKERR ufe.postProc("out_0.h5m");
448
449 int nb_levels = 5; // default number of refinement levels
450 // get number of refinement levels form command line
451 CHKERR PetscOptionsGetInt(PETSC_NULL, PETSC_NULL, "-nb_levels", &nb_levels,
452 PETSC_NULL);
453
454 // refine mesh, solve problem and do it again until number of refinement
455 // levels are exceeded.
456 for (int ll = 1; ll != nb_levels; ll++) {
457 const int nb_levels = ll;
458 CHKERR ufe.squashBits();
459 CHKERR ufe.refineMesh(ufe, nb_levels, order);
460 ref_level = BitRefLevel().set(nb_levels);
461 bc_flux_map.clear();
462 CHKERR ufe.addBoundaryElements(ref_level);
463 CHKERR ufe.buildProblem(ref_level);
467 CHKERR ufe.evaluateError();
469 CHKERR ufe.postProc(
470 static_cast<std::ostringstream &>(std::ostringstream().seekp(0)
471 << "out_" << nb_levels << ".h5m")
472 .str());
473 }
474 }
476
478
479 return 0;
480}
const std::string default_options
std::string param_file
Mix implementation of transport element.
int main()
Definition: adol-c_atom.cpp:46
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
#define BITREFLEVEL_SIZE
max number of refinements
Definition: definitions.h:219
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ TEMPERATURESET
Definition: definitions.h:155
@ HEATFLUXSET
Definition: definitions.h:156
@ NODESET
Definition: definitions.h:146
@ SIDESET
Definition: definitions.h:147
@ MAT_THERMALSET
block name is "MAT_THERMAL"
Definition: definitions.h:161
@ BLOCKSET
Definition: definitions.h:148
@ MOFEM_INVALID_DATA
Definition: definitions.h:36
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
virtual const RefEntity_multiIndex * get_ref_ents() const =0
Get the ref ents object.
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 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.
#define _IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet in a moFEM field.
#define _IT_CUBITMESHSETS_FOR_LOOP_(MESHSET_MANAGER, IT)
Iterator that loops over all the Cubit MeshSets in a moFEM field.
static char help[]
map< int, BcFluxData > BcFluxMap
auto bit
set bit
char mesh_file_name[255]
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
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 PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
MoFEMErrorCode postProc(const string out_file)
Post process results.
MoFEMErrorCode destroyMatrices()
destroy matrices
MoFEMErrorCode evaluateError()
Calculate error on elements.
MoFEMErrorCode solveLinearProblem()
solve problem
MoFEMErrorCode createMatrices()
create matrices
MoFEMErrorCode calculateResidual()
calculate residual
MoFEMErrorCode addFiniteElements(const std::string &fluxes_name, const std::string &values_name)
add finite elements
map< double, EntityHandle > errorMap
MoFEMErrorCode addFields(const std::string &values, const std::string &fluxes, const int order)
Add fields to database.
std::map< int, BlockData > setOfBlocks
maps block set id with appropriate BlockData
MoFEMErrorCode buildProblem(BitRefLevel &ref_level)
Build problem.
Managing BitRefLevels.
virtual moab::Interface & get_moab()=0
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.
Definition of the heat flux bc data structure.
Definition: BCData.hpp:422
Thermal material data structure.
Mesh refinement interface.
MoFEMErrorCode refineTets(const EntityHandle meshset, const BitRefLevel &bit, int verb=QUIET, const bool debug=false)
refine TET in the meshset
MoFEMErrorCode addVerticesInTheMiddleOfEdges(const EntityHandle meshset, const BitRefLevel &bit, const bool recursive=false, int verb=QUIET, EntityHandle start_v=0)
make vertices in the middle of edges in meshset and add them to refinement levels defined by bit
Interface for managing meshsets containing materials and boundary conditions.
MoFEMErrorCode setMeshsetFromFile(const string file_name, const bool clean_file_options=true)
add blocksets reading config file
Struct keeps handle to refined handle.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Application of mix transport data structure.
MoFEMErrorCode updateMeshsetsFieldsAndElements(const int nb_levels)
update meshsets with new entities after mesh refinement
MyTransport(MoFEM::Interface &m_field, BcFluxMap &bc_flux_map)
MoFEMErrorCode getBcOnValues(const EntityHandle ent, const double x, const double y, const double z, double &value)
natural (Dirihlet) boundary conditions (set values)
MoFEMErrorCode addBoundaryElements(BitRefLevel &ref_level)
set-up boundary conditions
MoFEMErrorCode getBcOnFluxes(const EntityHandle ent, const double x, const double y, const double z, double &flux)
essential (Neumann) boundary condition (set fluxes)
MoFEMErrorCode squashBits()
Squash bits of entities.
EntityHandle lastEnt
MoFEMErrorCode refineMesh(MixTransportElement &ufe, const int nb_levels, const int order)
Refine mesh.
MoFEMErrorCode getSource(EntityHandle ent, const double x, const double y, const double z, double &flux)
set source term
BcFluxMap & bcFluxMap