v0.6.20
Functions | Variables
thermal_steady.cpp File Reference

Example of steady thermal analysis. More...

#include <BasicFiniteElements.hpp>

Go to the source code of this file.

Functions

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

Variables

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

Detailed Description

Example of steady thermal analysis.

TODO:

Todo:
Make it work in distributed meshes with multigird solver. At the moment it is not working efficient as can.

Definition in file thermal_steady.cpp.

Function Documentation

◆ main()

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

Definition at line 34 of file thermal_steady.cpp.

34  {
35 
36  MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
37 
38  try {
39 
40  moab::Core mb_instance;
41  moab::Interface& moab = mb_instance;
42  int rank;
43  MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
44 
45  PetscBool flg = PETSC_TRUE;
46  char mesh_file_name[255];
47  CHKERR PetscOptionsGetString(PETSC_NULL, PETSC_NULL, "-my_file",
48  mesh_file_name, 255, &flg);
49  if (flg != PETSC_TRUE) {
50  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_FOUND,
51  "*** ERROR -my_file (MESH FILE NEEDED)");
52  }
53 
54  ParallelComm* pcomm = ParallelComm::get_pcomm(&moab,MYPCOMM_INDEX);
55  if(pcomm == NULL) pcomm = new ParallelComm(&moab,PETSC_COMM_WORLD);
56 
57  const char *option;
58  option = ""; //"PARALLEL=BCAST;";//;DEBUG_IO";
59  CHKERR moab.load_file(mesh_file_name, 0, option);
60 
61  // Create MoFEM (Joseph) database
62  MoFEM::Core core(moab);
63  MoFEM::Interface &m_field = core;
64 
65  //set entities bit level
66  BitRefLevel bit_level0;
67  bit_level0.set(0);
68  EntityHandle meshset_level0;
69  CHKERR moab.create_meshset(MESHSET_SET, meshset_level0);
70  CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevelByDim(0, 3,
71  bit_level0);
72 
73  // Fields
74  CHKERR m_field.add_field("TEMP", H1, AINSWORTH_LEGENDRE_BASE, 1);
75 
76  // Problem
77  CHKERR m_field.add_problem("THERMAL_PROBLEM");
78 
79  // set refinement level for problem
80  CHKERR m_field.modify_problem_ref_level_add_bit("THERMAL_PROBLEM",
81  bit_level0);
82 
83  // meshset consisting all entities in mesh
84  EntityHandle root_set = moab.get_root_set();
85  // add entities to field
86  CHKERR m_field.add_ents_to_field_by_type(root_set, MBTET, "TEMP");
87 
88  // set app. order
89  // see Hierarchic Finite Element Bases on Unstructured Tetrahedral Meshes
90  // (Mark Ainsworth & Joe Coyle)
91  PetscInt order;
92  CHKERR PetscOptionsGetInt(PETSC_NULL, PETSC_NULL, "-my_order", &order, &flg);
93  if (flg != PETSC_TRUE) {
94  order = 2;
95  }
96 
97  CHKERR m_field.set_field_order(root_set, MBTET, "TEMP", order);
98  CHKERR m_field.set_field_order(root_set, MBTRI, "TEMP", order);
99  CHKERR m_field.set_field_order(root_set, MBEDGE, "TEMP", order);
100  CHKERR m_field.set_field_order(root_set, MBVERTEX, "TEMP", 1);
101 
102  CHKERR m_field.add_field("MESH_NODE_POSITIONS", H1, AINSWORTH_LEGENDRE_BASE,
103  3);
104  CHKERR m_field.add_ents_to_field_by_type(root_set, MBTET,
105  "MESH_NODE_POSITIONS");
106  CHKERR m_field.set_field_order(0, MBTET, "MESH_NODE_POSITIONS", 2);
107  CHKERR m_field.set_field_order(0, MBTRI, "MESH_NODE_POSITIONS", 2);
108  CHKERR m_field.set_field_order(0, MBEDGE, "MESH_NODE_POSITIONS", 2);
109  CHKERR m_field.set_field_order(0, MBVERTEX, "MESH_NODE_POSITIONS", 1);
110 
111  ThermalElement thermal_elements(m_field);
112  CHKERR thermal_elements.addThermalElements("TEMP");
113  CHKERR thermal_elements.addThermalFluxElement("TEMP");
114  CHKERR thermal_elements.addThermalConvectionElement("TEMP");
115 
116  CHKERR m_field.modify_problem_add_finite_element("THERMAL_PROBLEM",
117  "THERMAL_FE");
118  CHKERR m_field.modify_problem_add_finite_element("THERMAL_PROBLEM",
119  "THERMAL_FLUX_FE");
120  CHKERR m_field.modify_problem_add_finite_element("THERMAL_PROBLEM",
121  "THERMAL_CONVECTION_FE");
122 
123  /****/
124  // build database
125  // build field
126  CHKERR m_field.build_fields();
127  // build finite elemnts
128  CHKERR m_field.build_finite_elements();
129  // build adjacencies
130  CHKERR m_field.build_adjacencies(bit_level0);
131 
132  ProblemsManager *prb_mng_ptr;
133  CHKERR m_field.getInterface(prb_mng_ptr);
134  // build problem
135  CHKERR prb_mng_ptr->buildProblem("THERMAL_PROBLEM", true);
136 
137  Projection10NodeCoordsOnField ent_method_material(m_field,
138  "MESH_NODE_POSITIONS");
139  CHKERR m_field.loop_dofs("MESH_NODE_POSITIONS", ent_method_material);
140 
141  /****/
142  //mesh partitioning
143  //partition
144  CHKERR prb_mng_ptr->partitionProblem("THERMAL_PROBLEM");
145  CHKERR prb_mng_ptr->partitionFiniteElements("THERMAL_PROBLEM");
146  // what are ghost nodes, see Petsc Manual
147  CHKERR prb_mng_ptr->partitionGhostDofs("THERMAL_PROBLEM");
148 
149  Vec F;
150  CHKERR m_field.getInterface<VecManager>()->vecCreateGhost("THERMAL_PROBLEM",
151  ROW, &F);
152  Vec T;
153  CHKERR VecDuplicate(F,&T);
154  Mat A;
155  CHKERR m_field.MatCreateMPIAIJWithArrays("THERMAL_PROBLEM",&A);
156 
157  DirichletTemperatureBc my_dirichlet_bc(m_field, "TEMP", A, T, F);
158  CHKERR thermal_elements.setThermalFiniteElementRhsOperators("TEMP", F);
159  CHKERR thermal_elements.setThermalFiniteElementLhsOperators("TEMP", A);
160  CHKERR thermal_elements.setThermalFluxFiniteElementRhsOperators("TEMP", F);
161  CHKERR thermal_elements.setThermalConvectionFiniteElementLhsOperators("TEMP",
162  A);
163  CHKERR thermal_elements.setThermalConvectionFiniteElementRhsOperators("TEMP",
164  F);
165 
166  CHKERR VecZeroEntries(T);
167  CHKERR VecGhostUpdateBegin(T, INSERT_VALUES, SCATTER_FORWARD);
168  CHKERR VecGhostUpdateEnd(T, INSERT_VALUES, SCATTER_FORWARD);
169  CHKERR VecZeroEntries(F);
170  CHKERR VecGhostUpdateBegin(F, INSERT_VALUES, SCATTER_FORWARD);
171  CHKERR VecGhostUpdateEnd(F, INSERT_VALUES, SCATTER_FORWARD);
172  CHKERR MatZeroEntries(A);
173 
174  //preproc
175  CHKERR m_field.problem_basic_method_preProcess("THERMAL_PROBLEM",
176  my_dirichlet_bc);
177  CHKERR m_field.getInterface<VecManager>()->setGlobalGhostVector(
178  "THERMAL_PROBLEM", ROW, T, INSERT_VALUES, SCATTER_REVERSE);
179 
180  CHKERR m_field.loop_finite_elements("THERMAL_PROBLEM", "THERMAL_FE",
181  thermal_elements.getLoopFeRhs());
182  CHKERR m_field.loop_finite_elements("THERMAL_PROBLEM", "THERMAL_FE",
183  thermal_elements.getLoopFeLhs());
184  CHKERR m_field.loop_finite_elements("THERMAL_PROBLEM", "THERMAL_FLUX_FE",
185  thermal_elements.getLoopFeFlux());
187  "THERMAL_PROBLEM", "THERMAL_CONVECTION_FE",
188  thermal_elements.getLoopFeConvectionRhs());
190  "THERMAL_PROBLEM", "THERMAL_CONVECTION_FE",
191  thermal_elements.getLoopFeConvectionLhs());
192 
193  //postproc
194  CHKERR m_field.problem_basic_method_postProcess("THERMAL_PROBLEM",
195  my_dirichlet_bc);
196 
197  CHKERR VecGhostUpdateBegin(F, ADD_VALUES, SCATTER_REVERSE);
198  CHKERR VecGhostUpdateEnd(F, ADD_VALUES, SCATTER_REVERSE);
199  CHKERR VecAssemblyBegin(F);
200  CHKERR VecAssemblyEnd(F);
201  CHKERR MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
202  CHKERR MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
203 
204  CHKERR VecScale(F, -1);
205 
206  // Solver
207  KSP solver;
208  CHKERR KSPCreate(PETSC_COMM_WORLD, &solver);
209  CHKERR KSPSetOperators(solver, A, A);
210  CHKERR KSPSetFromOptions(solver);
211  CHKERR KSPSetUp(solver);
212 
213  CHKERR KSPSolve(solver, F, T);
214  CHKERR VecGhostUpdateBegin(T, INSERT_VALUES, SCATTER_FORWARD);
215  CHKERR VecGhostUpdateEnd(T, INSERT_VALUES, SCATTER_FORWARD);
216 
217  CHKERR m_field.problem_basic_method_preProcess("THERMAL_PROBLEM",
218  my_dirichlet_bc);
219 
220  //Save data on mesh
221  CHKERR m_field.getInterface<VecManager>()->setGlobalGhostVector(
222  "THERMAL_PROBLEM", ROW, T, INSERT_VALUES, SCATTER_REVERSE);
223  //CHKERR VecView(F,PETSC_VIEWER_STDOUT_WORLD);
224 
225  //Range ref_edges;
226  //CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(bit_level0,BitRefLevel().set(),MBEDGE,ref_edges);
227  //CHKERR moab.list_entities(ref_edges);
228  //m_field.list_dofs_by_field_name("TEMP");
229 
230  if(pcomm->rank()==0) {
231  CHKERR moab.write_file("solution.h5m");
232  }
233 
234  /*EntityHandle fe_meshset = m_field.get_finite_element_meshset("THERMAL_FE");
235  Range tets;
236  CHKERR moab.get_entities_by_type(fe_meshset,MBTET,tets,true);
237  Range tets_edges;
238  CHKERR moab.get_adjacencies(tets,1,false,tets_edges,moab::Interface::UNION);
239  EntityHandle edges_meshset;
240  CHKERR moab.create_meshset(MESHSET_SET,edges_meshset);
241  CHKERR moab.add_entities(edges_meshset,tets);
242  CHKERR moab.add_entities(edges_meshset,tets_edges);
243  CHKERR moab.convert_entities(edges_meshset,true,false,false); */
244 
245  ProjectionFieldOn10NodeTet ent_method_on_10nodeTet(m_field, "TEMP", true,
246  false, "TEMP");
247  CHKERR m_field.loop_dofs("TEMP", ent_method_on_10nodeTet);
248  ent_method_on_10nodeTet.setNodes = false;
249  CHKERR m_field.loop_dofs("TEMP", ent_method_on_10nodeTet);
250 
251  if (pcomm->rank() == 0) {
252  EntityHandle out_meshset;
253  CHKERR moab.create_meshset(MESHSET_SET, out_meshset);
255  "THERMAL_PROBLEM", "THERMAL_FE", out_meshset);
256  CHKERR moab.write_file("out.vtk", "VTK", "", &out_meshset, 1);
257  CHKERR moab.delete_entities(&out_meshset, 1);
258  }
259 
260  CHKERR MatDestroy(&A);
261  CHKERR VecDestroy(&F);
262  CHKERR VecDestroy(&T);
263  CHKERR KSPDestroy(&solver);
264 
265  }
266  CATCH_ERRORS;
267 
268  return MoFEM::Core::Finalize();
269 }
Problem manager is used to build and partition problems .
virtual MoFEMErrorCode add_problem(const std::string &name, enum MoFEMTypes bh=MF_EXCL, int verb=-1)=0
Add 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 ...
Core (interface) class.
Definition: Core.hpp:50
Projection of edge entities with one mid-node on hierarchical basis.
virtual MoFEMErrorCode build_finite_elements(int verb=-1)=0
Build finite elementsBuild finite element data structures. Have to be run before problem and adjacenc...
virtual MoFEMErrorCode loop_finite_elements(const std::string &problem_name, const std::string &fe_name, FEMethod &method, MoFEMTypes bh=MF_EXIST, int verb=-1)=0
Make a loop over finite elements.
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
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=-1)=0
Add field.
MoFEMErrorCode buildProblem(const std::string &name, const bool square_matrix, int verb=1)
build problem data structures
virtual MoFEMErrorCode MatCreateMPIAIJWithArrays(const std::string &name, Mat *Aij, int verb=-1)=0
create Mat (MPIAIJ) for problem (collective)
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Common.hpp:147
MoFEMErrorCode getInterface(const MOFEMuuid &uuid, IFACE *&iface) const
Get interface by uuid and return reference to pointer of interface.
MoFEMErrorCode partitionFiniteElements(const std::string &name, bool part_from_moab=false, int low_proc=-1, int hi_proc=-1, int verb=1)
partition finite elementsFunction which partition finite elements based on dofs partitioning. In addition it sets information about local row and cols dofs at given element on partition.
virtual MoFEMErrorCode add_ents_to_field_by_type(const Range &ents, const EntityType type, const std::string &name, int verb=-1)=0
Add entities to field meshsetThe lower dimension entities are added depending on the space type...
Managing BitRefLevels.
virtual MoFEMErrorCode set_field_order(const EntityHandle meshset, const EntityType type, const std::string &name, const ApproximationOrder order, int verb=-1)=0
Set order approximation of the entities in the field.
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:141
static char help[]
virtual MoFEMErrorCode modify_problem_ref_level_add_bit(const std::string &name_problem, const BitRefLevel &bit)=0
add ref level to problemif same finite element is solved using different level of refinements...
Vector manager is used to create vectors .
Definition: VecManager.hpp:35
#define CHKERR
Inline error check.
Definition: definitions.h:610
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:316
MoFEMErrorCode partitionGhostDofs(const std::string &name, int verb=1)
determine ghost nodes
virtual MoFEMErrorCode build_fields(int verb=-1)=0
virtual MoFEMErrorCode loop_dofs(const Problem *problem_ptr, const std::string &field_name, RowColData rc, EntMethod &method, int lower_rank, int upper_rank, int verb=-1)=0
Make a loop over entities.
virtual MoFEMErrorCode get_problem_finite_elements_entities(const std::string &name, const std::string &fe_name, const EntityHandle meshset)=0
add finite elements to the meshset
MoFEMErrorCode partitionProblem(const std::string &name, int verb=1)
partition problem dofs (collective)
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:462
continuous field
Definition: definitions.h:176
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=-1)=0
build adjacencies
virtual MoFEMErrorCode problem_basic_method_preProcess(const Problem *problem_ptr, BasicMethod &method, int verb=-1)=0
Set data for BasicMethod.
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:61
structure grouping operators and data used for thermal problemsIn order to assemble matrices and righ...
virtual MoFEMErrorCode problem_basic_method_postProcess(const Problem *problem_ptr, BasicMethod &method, int verb=-1)=0
Set data for BasicMethodThis function set data about problem, adjacencies and other multi-indices in ...

Variable Documentation

◆ help

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

Definition at line 32 of file thermal_steady.cpp.