v0.7.26
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 build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
MoFEMErrorCode buildProblem(const std::string &name, const bool square_matrix, int verb=VERBOSE)
build problem data structures
MoFEMErrorCode partitionGhostDofs(const std::string &name, int verb=VERBOSE)
determine ghost nodes
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. In addition it sets information about local row and cols dofs at given element on partition.
MoFEMErrorCode partitionProblem(const std::string &name, int verb=VERBOSE)
partition problem dofs (collective)
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 ...
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
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
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
virtual MoFEMErrorCode loop_finite_elements(const std::string &problem_name, const std::string &fe_name, FEMethod &method, MoFEMTypes bh=MF_EXIST, int verb=DEFAULT_VERBOSITY)=0
Make a loop over finite elements.
MoFEMErrorCode getInterface(const MOFEMuuid &uuid, IFACE *&iface) const
Get interface by uuid and return reference to pointer of interface.
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:142
static char help[]
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
Vector manager is used to create vectors .
Definition: VecManager.hpp:35
virtual MoFEMErrorCode problem_basic_method_preProcess(const Problem *problem_ptr, BasicMethod &method, int verb=DEFAULT_VERBOSITY)=0
Set data for BasicMethod.
#define CHKERR
Inline error check.
Definition: definitions.h:617
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
virtual MoFEMErrorCode MatCreateMPIAIJWithArrays(const std::string &name, Mat *Aij, int verb=DEFAULT_VERBOSITY)=0
create Mat (MPIAIJ) for problem (collective)
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:324
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.
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
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:469
continuous field
Definition: definitions.h:177
virtual MoFEMErrorCode problem_basic_method_postProcess(const Problem *problem_ptr, BasicMethod &method, int verb=DEFAULT_VERBOSITY)=0
Set data for BasicMethodThis function set data about problem, adjacencies and other multi-indices in ...
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
structure grouping operators and data used for thermal problemsIn order to assemble matrices and righ...

Variable Documentation

◆ help

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

Definition at line 32 of file thermal_steady.cpp.