v0.14.0
Loading...
Searching...
No Matches
gel_analysis.cpp
Go to the documentation of this file.
1/** \file gel_analysis.cpp
2 \brief Reads cubit file and solves problem with gel material
3 \ingroup gel
4
5 \todo Current version is limited only to one material. If one like to have
6 general problem for nonlinear elasticity should implement general time
7 dependent problem. If inertia terms need to be considered, this material
8 should be add to nonlinear dynamics problem.
9
10 \todo Internal history state variables need to be statically condensed. It
11 can be done by implementing static condensation on finite element level or
12 by implementing pre-conditioner.
13
14*/
15
16/* This file is part of MoFEM.
17 * MoFEM is free software: you can redistribute it and/or modify it under
18 * the terms of the GNU Lesser General Public License as published by the
19 * Free Software Foundation, either version 3 of the License, or (at your
20 * option) any later version.
21 *
22 * MoFEM is distributed in the hope that it will be useful, but WITHOUT
23 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
24 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
25 * License for more details.
26 *
27 * You should have received a copy of the GNU Lesser General Public
28 * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
29
30#include <MoFEM.hpp>
31using namespace MoFEM;
32
33#include <PostProcOnRefMesh.hpp>
35
36#include <boost/numeric/ublas/vector_proxy.hpp>
37#include <boost/numeric/ublas/matrix.hpp>
38#include <boost/numeric/ublas/matrix_proxy.hpp>
39#include <boost/numeric/ublas/vector.hpp>
40#include <adolc/adolc.h>
41#include <Gels.hpp>
42#include <UserGelModel.hpp>
43using namespace GelModule;
44
45#include <boost/program_options.hpp>
46using namespace std;
47namespace po = boost::program_options;
48
50#include <DirichletBC.hpp>
51#include <SurfacePressure.hpp>
52#include <NodalForce.hpp>
53#include <EdgeForce.hpp>
54
55// Elements for applying fluxes, convection or radiation
56// are used to apply solvent concentration.
57#include <ThermalElement.hpp>
58
59//See file users_modules/elasticity/TimeForceScale.hpp
60#include <TimeForceScale.hpp>
61
62#include <boost/iostreams/tee.hpp>
63#include <boost/iostreams/stream.hpp>
64#include <fstream>
65#include <iostream>
66namespace bio = boost::iostreams;
67using bio::tee_device;
68using bio::stream;
69
70
71
72static char help[] = "...\n\n";
73
74struct BlockData {
75 int oRder;
77 oRder(-1) {
78 }
79};
80
81int main(int argc, char *argv[]) {
82
83 MoFEM::Core::Initialize(&argc,&argv,(char *)0,help);
84
85 try {
86
87 moab::Core mb_instance;
88 moab::Interface& moab = mb_instance;
89
90 PetscBool flg_gel_config,flg_file;
91 char mesh_file_name[255];
92 char gel_config_file[255];
93 int order = 2;
94 PetscBool is_partitioned = PETSC_FALSE;
95 PetscBool is_atom_test = PETSC_FALSE;
96
97 {
98
99 ierr = PetscOptionsBegin(PETSC_COMM_WORLD,"","Elastic Config","none"); CHKERRQ(ierr);
100 ierr = PetscOptionsString(
101 "-my_file",
102 "mesh file name","",
103 "mesh.h5m",
105 255,
106 &flg_file
107 ); CHKERRQ(ierr);
108 ierr = PetscOptionsInt(
109 "-my_order",
110 "default approximation order",
111 "",
112 order,
113 &order,
114 PETSC_NULL
115 ); CHKERRQ(ierr);
116
117 ierr = PetscOptionsBool(
118 "-my_is_partitioned",
119 "set if mesh is partitioned (this result that each process keeps only part of the mes","",
120 PETSC_FALSE,&is_partitioned,PETSC_NULL
121 ); CHKERRQ(ierr);
122 ierr = PetscOptionsString(
123 "-my_gel_config",
124 "gel configuration file name","",
125 "gel_config.in",gel_config_file,255,&flg_gel_config
126 ); CHKERRQ(ierr);
127
128 ierr = PetscOptionsBool(
129 "-my_is_atom_test",
130 "is used with testing, exit with error when diverged","",
131 PETSC_FALSE,&is_atom_test,PETSC_NULL
132 ); CHKERRQ(ierr);
133
134 ierr = PetscOptionsEnd(); CHKERRQ(ierr);
135
136 //Reade parameters from line command
137 if(flg_file != PETSC_TRUE) {
138 SETERRQ(PETSC_COMM_SELF,1,"*** ERROR -my_file (MESH FILE NEEDED)");
139 }
140
141 if(is_partitioned == PETSC_TRUE) {
142 //Read mesh to MOAB
143 const char *option;
144 option =
145 "PARALLEL=BCAST_DELETE;"
146 "PARALLEL_RESOLVE_SHARED_ENTS;"
147 "PARTITION=PARALLEL_PARTITION;";
148 rval = moab.load_file(mesh_file_name, 0, option); CHKERRQ_MOAB(rval);
149 } else {
150 const char *option;
151 option = "";
152 rval = moab.load_file(mesh_file_name, 0, option); CHKERRQ_MOAB(rval);
153 }
154 ParallelComm* pcomm = ParallelComm::get_pcomm(&moab,MYPCOMM_INDEX);
155 if(pcomm == NULL) pcomm = new ParallelComm(&moab,PETSC_COMM_WORLD);
156 }
157
158 MoFEM::Core core(moab);
159 MoFEM::Interface& m_field = core;
160
161 // Seed all mesh entities to MoFEM database, those entities can be potentially used as finite elements
162 // or as entities which carry some approximation field.
163 BitRefLevel bit_level0;
164 bit_level0.set(0);
165 ierr = m_field.seed_ref_level_3D(0,bit_level0); CHKERRQ(ierr);
166
167 // Define fields and finite elements
168 Gel gel(m_field);
169 map<int,ThermalElement::FluxData> set_of_solvent_fluxes;
170 {
171
172 // Set approximation fields
173 {
174
175 // Add fields
176 bool check_if_spatial_field_exist = m_field.check_field("SPATIAL_POSITION");
177 ierr = m_field.add_field("SPATIAL_POSITION",H1,AINSWORTH_LEGENDRE_BASE,3,MB_TAG_SPARSE,MF_ZERO); CHKERRQ(ierr);
178 ierr = m_field.add_field("CHEMICAL_LOAD",H1,AINSWORTH_LEGENDRE_BASE,1,MB_TAG_SPARSE,MF_ZERO); CHKERRQ(ierr);
179 ierr = m_field.add_field("HAT_EPS",L2,AINSWORTH_LEGENDRE_BASE,6,MB_TAG_SPARSE,MF_ZERO); CHKERRQ(ierr);
180 ierr = m_field.add_field("SPATIAL_POSITION_DOT",H1,AINSWORTH_LEGENDRE_BASE,3,MB_TAG_SPARSE,MF_ZERO); CHKERRQ(ierr);
181 ierr = m_field.add_field("CHEMICAL_LOAD_DOT",H1,AINSWORTH_LEGENDRE_BASE,1,MB_TAG_SPARSE,MF_ZERO); CHKERRQ(ierr);
182 ierr = m_field.add_field("HAT_EPS_DOT",L2,AINSWORTH_LEGENDRE_BASE,6,MB_TAG_SPARSE,MF_ZERO); CHKERRQ(ierr);
183
184 //Add field H1 space rank 3 to approximate geometry using hierarchical basis
185 //For 10 node tetrahedral, before use, geometry is projected on that field (see below)
186 ierr = m_field.add_field("MESH_NODE_POSITIONS",H1,AINSWORTH_LEGENDRE_BASE,3,MB_TAG_SPARSE,MF_ZERO); CHKERRQ(ierr);
187
188 //meshset consisting all entities in mesh
189 EntityHandle root_set = moab.get_root_set();
190 //Add entities to field (root_mesh, i.e. on all mesh entities fields are approx.)
191 ierr = m_field.add_ents_to_field_by_type(root_set,MBTET,"SPATIAL_POSITION"); CHKERRQ(ierr);
192 ierr = m_field.add_ents_to_field_by_type(root_set,MBTET,"CHEMICAL_LOAD"); CHKERRQ(ierr);
193 ierr = m_field.add_ents_to_field_by_type(root_set,MBTET,"HAT_EPS"); CHKERRQ(ierr);
194 ierr = m_field.add_ents_to_field_by_type(root_set,MBTET,"SPATIAL_POSITION_DOT"); CHKERRQ(ierr);
195 ierr = m_field.add_ents_to_field_by_type(root_set,MBTET,"CHEMICAL_LOAD_DOT"); CHKERRQ(ierr);
196 ierr = m_field.add_ents_to_field_by_type(root_set,MBTET,"HAT_EPS_DOT"); CHKERRQ(ierr);
197
198 // Set approximation order. Solvent concentration has one order less than
199 // order of of spatial position field. Tests need to be maid if that is
200 // good enough to have stability for this type of problem. If not bubble
201 // functions could be easily added by increasing approximate order for
202 // volume.
203
204 ierr = m_field.set_field_order(root_set,MBTET,"SPATIAL_POSITION",order,2); CHKERRQ(ierr);
205 ierr = m_field.set_field_order(root_set,MBTRI,"SPATIAL_POSITION",order,2); CHKERRQ(ierr);
206 ierr = m_field.set_field_order(root_set,MBEDGE,"SPATIAL_POSITION",order,2); CHKERRQ(ierr);
207 ierr = m_field.set_field_order(root_set,MBVERTEX,"SPATIAL_POSITION",1); CHKERRQ(ierr);
208
209 ierr = m_field.set_field_order(root_set,MBTET,"SPATIAL_POSITION_DOT",order); CHKERRQ(ierr);
210 ierr = m_field.set_field_order(root_set,MBTRI,"SPATIAL_POSITION_DOT",order); CHKERRQ(ierr);
211 ierr = m_field.set_field_order(root_set,MBEDGE,"SPATIAL_POSITION_DOT",order); CHKERRQ(ierr);
212 ierr = m_field.set_field_order(root_set,MBVERTEX,"SPATIAL_POSITION_DOT",1); CHKERRQ(ierr);
213
214 ierr = m_field.set_field_order(root_set,MBTET,"CHEMICAL_LOAD",order-1); CHKERRQ(ierr);
215 ierr = m_field.set_field_order(root_set,MBTRI,"CHEMICAL_LOAD",order-1); CHKERRQ(ierr);
216 ierr = m_field.set_field_order(root_set,MBEDGE,"CHEMICAL_LOAD",order-1); CHKERRQ(ierr);
217 ierr = m_field.set_field_order(root_set,MBVERTEX,"CHEMICAL_LOAD",1); CHKERRQ(ierr);
218
219 ierr = m_field.set_field_order(root_set,MBTET,"CHEMICAL_LOAD_DOT",order-1); CHKERRQ(ierr);
220 ierr = m_field.set_field_order(root_set,MBTRI,"CHEMICAL_LOAD_DOT",order-1); CHKERRQ(ierr);
221 ierr = m_field.set_field_order(root_set,MBEDGE,"CHEMICAL_LOAD_DOT",order-1); CHKERRQ(ierr);
222 ierr = m_field.set_field_order(root_set,MBVERTEX,"CHEMICAL_LOAD_DOT",1); CHKERRQ(ierr);
223
224 ierr = m_field.set_field_order(root_set,MBTET,"HAT_EPS",order-1); CHKERRQ(ierr);
225 ierr = m_field.set_field_order(root_set,MBTET,"HAT_EPS_DOT",order-1); CHKERRQ(ierr);
226
227 //gemetry approximation is set to 2nd oreder
228 ierr = m_field.add_ents_to_field_by_type(root_set,MBTET,"MESH_NODE_POSITIONS"); CHKERRQ(ierr);
229 ierr = m_field.set_field_order(0,MBTET,"MESH_NODE_POSITIONS",2); CHKERRQ(ierr);
230 ierr = m_field.set_field_order(0,MBTRI,"MESH_NODE_POSITIONS",2); CHKERRQ(ierr);
231 ierr = m_field.set_field_order(0,MBEDGE,"MESH_NODE_POSITIONS",2); CHKERRQ(ierr);
232 ierr = m_field.set_field_order(0,MBVERTEX,"MESH_NODE_POSITIONS",1); CHKERRQ(ierr);
233
234 try {
235
236 map<int,BlockData> block_data;
237 map<int,Gel::BlockMaterialData> &material_data = gel.blockMaterialData;
238
239 ifstream ini_file(gel_config_file);
240 po::variables_map vm;
241 po::options_description config_file_options;
242
243 int def_block = 0;
244 Tag th_block_id;
245 rval = moab.tag_get_handle(
246 "BLOCK_ID",1,MB_TYPE_INTEGER,th_block_id,MB_TAG_CREAT|MB_TAG_SPARSE,&def_block
247 ); CHKERRQ_MOAB(rval);
248
250
251 const string &name = it->getName();
252 if(name.compare(0,3,"GEL")!=0) continue;
253
254 int block_id = it->getMeshsetId();
255
256 ostringstream str_order;
257 str_order << "block_" << block_id << ".oRder";
258 config_file_options.add_options()
259 (str_order.str().c_str(),po::value<int>(
260 &block_data[block_id].oRder
261 )->default_value(order));
262 Range block_ents;
263 rval = moab.get_entities_by_handle(
264 it->getMeshset(),block_ents,true
265 ); CHKERRQ_MOAB(rval);
266 material_data[block_id].tEts = block_ents.subset_by_type(MBTET);
267 std::vector<int> block_id_vec(material_data[block_id].tEts.size(),it->getMeshsetId());
268 // cerr << "AAA " << block_id_vec.size() << " " << material_data[block_id].tEts.size() << endl;
269 rval = moab.tag_set_data(
270 th_block_id,material_data[block_id].tEts,&*block_id_vec.begin()
271 );
272
273 ostringstream str_g_alpha;
274 str_g_alpha << "block_" << block_id << ".gAlpha";
275 config_file_options.add_options()
276 (str_g_alpha.str().c_str(),po::value<double>(&material_data[block_id].gAlpha)->default_value(1));
277 ostringstream str_v_alpha;
278 str_v_alpha << "block_" << block_id << ".vAlpha";
279 config_file_options.add_options()
280 (str_v_alpha.str().c_str(),po::value<double>(&material_data[block_id].vAlpha)->default_value(0));
281 ostringstream str_g_beta;
282 str_g_beta << "block_" << block_id << ".gBeta";
283 config_file_options.add_options()
284 (str_g_beta.str().c_str(),po::value<double>(&material_data[block_id].gBeta)->default_value(1));
285 ostringstream str_v_beta;
286 str_v_beta << "block_" << block_id << ".vBeta";
287 config_file_options.add_options()
288 (str_v_beta.str().c_str(),po::value<double>(&material_data[block_id].vBeta)->default_value(0));
289 ostringstream str_g_beta_hat;
290 str_g_beta_hat << "block_" << block_id << ".gBetaHat";
291 config_file_options.add_options()
292 (str_g_beta_hat.str().c_str(),po::value<double>(&material_data[block_id].gBetaHat)->default_value(1));
293 ostringstream str_v_beta_hat;
294 str_v_beta_hat << "block_" << block_id << ".vBetaHat";
295 config_file_options.add_options()
296 (str_v_beta_hat.str().c_str(),po::value<double>(&material_data[block_id].vBetaHat)->default_value(0));
297 ostringstream str_omega;
298 str_omega << "block_" << block_id << ".oMega";
299 config_file_options.add_options()
300 (str_omega.str().c_str(),po::value<double>(&material_data[block_id].oMega)->default_value(1));
301 ostringstream str_viscosity;
302 str_viscosity << "block_" << block_id << ".vIscosity";
303 config_file_options.add_options()
304 (str_viscosity.str().c_str(),po::value<double>(&material_data[block_id].vIscosity)->default_value(1));
305 ostringstream str_permeability;
306 str_permeability << "block_" << block_id << ".pErmeability";
307 config_file_options.add_options()
308 (str_permeability.str().c_str(),po::value<double>(&material_data[block_id].pErmeability)->default_value(1));
309 ostringstream str_mu0;
310 str_mu0 << "block_" << block_id << ".mU0";
311 config_file_options.add_options()
312 (str_mu0.str().c_str(),po::value<double>(&material_data[block_id].mU0)->default_value(1));
313
314 }
315
316 po::parsed_options parsed = parse_config_file(ini_file,config_file_options,true);
317 store(parsed,vm);
318 po::notify(vm);
320
321 const string &name = it->getName();
322 if(name.compare(0,3,"GEL")!=0) continue;
323
324 if(block_data[it->getMeshsetId()].oRder == -1) continue;
325 if(block_data[it->getMeshsetId()].oRder == order) continue;
326 PetscPrintf(PETSC_COMM_WORLD,"Set block %d order to %d\n",it->getMeshsetId(),block_data[it->getMeshsetId()].oRder);
327 Range block_ents;
328 rval = moab.get_entities_by_handle(it->getMeshset(),block_ents,true); CHKERRQ_MOAB(rval);
329 Range ents_to_set_order;
330 rval = moab.get_adjacencies(block_ents,3,false,ents_to_set_order,moab::Interface::UNION); CHKERRQ_MOAB(rval);
331 ents_to_set_order = ents_to_set_order.subset_by_type(MBTET);
332 rval = moab.get_adjacencies(block_ents,2,false,ents_to_set_order,moab::Interface::UNION); CHKERRQ_MOAB(rval);
333 rval = moab.get_adjacencies(block_ents,1,false,ents_to_set_order,moab::Interface::UNION); CHKERRQ_MOAB(rval);
334 ierr = m_field.synchronise_entities(ents_to_set_order); CHKERRQ(ierr);
335 ierr = m_field.set_field_order(ents_to_set_order,"SPATIAL_POSITION",block_data[it->getMeshsetId()].oRder); CHKERRQ(ierr);
336 ierr = m_field.set_field_order(ents_to_set_order,"CHEMICAL_LOAD",block_data[it->getMeshsetId()].oRder-1); CHKERRQ(ierr);
337 ierr = m_field.set_field_order(ents_to_set_order,"HAT_EPS",block_data[it->getMeshsetId()].oRder-1); CHKERRQ(ierr);
338 ierr = m_field.set_field_order(ents_to_set_order,"SPATIAL_POSITION_DOT",block_data[it->getMeshsetId()].oRder); CHKERRQ(ierr);
339 ierr = m_field.set_field_order(ents_to_set_order,"CHEMICAL_LOAD_DOT",block_data[it->getMeshsetId()].oRder-1); CHKERRQ(ierr);
340 ierr = m_field.set_field_order(ents_to_set_order,"HAT_EPS_DOT",block_data[it->getMeshsetId()].oRder-1); CHKERRQ(ierr);
341 }
342 vector<string> additional_parameters;
343 additional_parameters = collect_unrecognized(parsed.options,po::include_positional);
344 for(vector<string>::iterator vit = additional_parameters.begin();
345 vit!=additional_parameters.end();vit++) {
346 ierr = PetscPrintf(PETSC_COMM_WORLD,"** WARNING Unrecognised option %s\n",vit->c_str()); CHKERRQ(ierr);
347 }
348 } catch (const std::exception& ex) {
349 ostringstream ss;
350 ss << ex.what() << endl;
351 SETERRQ(PETSC_COMM_SELF,MOFEM_STD_EXCEPTION_THROW,ss.str().c_str());
352 }
353
354 for(
355 map<int,Gel::BlockMaterialData>::iterator mit = gel.blockMaterialData.begin();
356 mit!=gel.blockMaterialData.end();mit++
357 ) {
358 ierr = PetscPrintf(PETSC_COMM_WORLD,"** Block Id %d\n",mit->first); CHKERRQ(ierr);
359 PetscPrintf(PETSC_COMM_WORLD,"vAlpha = %5.4g\n",mit->second.vAlpha);
360 PetscPrintf(PETSC_COMM_WORLD,"gAlpha = %5.4g\n",mit->second.gAlpha);
361 PetscPrintf(PETSC_COMM_WORLD,"vBeta = %5.4g\n",mit->second.vBeta);
362 PetscPrintf(PETSC_COMM_WORLD,"gBeta = %5.4g\n",mit->second.gBeta);
363 PetscPrintf(PETSC_COMM_WORLD,"vBetaHat = %5.4g\n",mit->second.vBetaHat);
364 PetscPrintf(PETSC_COMM_WORLD,"gBetaHat = %5.4g\n",mit->second.gBetaHat);
365 PetscPrintf(PETSC_COMM_WORLD,"pErmeability = %5.4g\n",mit->second.pErmeability);
366 PetscPrintf(PETSC_COMM_WORLD,"vIscosity = %5.4g\n",mit->second.vIscosity);
367 PetscPrintf(PETSC_COMM_WORLD,"oMega = %5.4g\n",mit->second.oMega);
368 PetscPrintf(PETSC_COMM_WORLD,"mU0 = %5.4g\n",mit->second.mU0);
369 }
370
371 ierr = m_field.build_fields(); CHKERRQ(ierr);
372
373 // Sett geometry approximation and initial spatial positions
374 // 10 node tets
375 if (!check_if_spatial_field_exist) {
376 Projection10NodeCoordsOnField ent_method_material(m_field, "MESH_NODE_POSITIONS");
377 ierr = m_field.loop_dofs("MESH_NODE_POSITIONS", ent_method_material); CHKERRQ(ierr);
378 Projection10NodeCoordsOnField ent_method_spatial(m_field, "SPATIAL_POSITION");
379 ierr = m_field.loop_dofs("SPATIAL_POSITION", ent_method_spatial); CHKERRQ(ierr);
380 // Set value to chemical load
381 map<int,Gel::BlockMaterialData> &material_data = gel.blockMaterialData;
382 for(
383 map<int,Gel::BlockMaterialData>::iterator mit = material_data.begin();
384 mit!=material_data.end();mit++
385 ) {
386 // FIXME mu0 at border in undefined, if two blackset have different mu0
387 Range vertices;
388 rval = moab.get_connectivity(mit->second.tEts,vertices,true); CHKERRQ_MOAB(rval);
389 ierr = m_field.set_field(mit->second.mU0,MBVERTEX,vertices,"CHEMICAL_LOAD"); CHKERRQ(ierr);
390 }
391 }
392
393 }
394
395 //Set finite elements. The primary element is GEL_FE in addition elements
396 //for applying tractions and fluxes are added.
397 {
398 ierr = m_field.add_finite_element("GEL_FE",MF_ZERO); CHKERRQ(ierr);
399 ierr = m_field.modify_finite_element_add_field_row("GEL_FE","SPATIAL_POSITION"); CHKERRQ(ierr);
400 ierr = m_field.modify_finite_element_add_field_col("GEL_FE","SPATIAL_POSITION"); CHKERRQ(ierr);
401 ierr = m_field.modify_finite_element_add_field_row("GEL_FE","CHEMICAL_LOAD"); CHKERRQ(ierr);
402 ierr = m_field.modify_finite_element_add_field_col("GEL_FE","CHEMICAL_LOAD"); CHKERRQ(ierr);
403 ierr = m_field.modify_finite_element_add_field_row("GEL_FE","HAT_EPS"); CHKERRQ(ierr);
404 ierr = m_field.modify_finite_element_add_field_col("GEL_FE","HAT_EPS"); CHKERRQ(ierr);
405 ierr = m_field.modify_finite_element_add_field_data("GEL_FE","SPATIAL_POSITION"); CHKERRQ(ierr);
406 ierr = m_field.modify_finite_element_add_field_data("GEL_FE","SPATIAL_POSITION_DOT"); CHKERRQ(ierr);
407 ierr = m_field.modify_finite_element_add_field_data("GEL_FE","CHEMICAL_LOAD"); CHKERRQ(ierr);
408 ierr = m_field.modify_finite_element_add_field_data("GEL_FE","CHEMICAL_LOAD_DOT"); CHKERRQ(ierr);
409 ierr = m_field.modify_finite_element_add_field_data("GEL_FE","HAT_EPS"); CHKERRQ(ierr);
410 ierr = m_field.modify_finite_element_add_field_data("GEL_FE","HAT_EPS_DOT"); CHKERRQ(ierr);
411 ierr = m_field.modify_finite_element_add_field_data("GEL_FE","MESH_NODE_POSITIONS"); CHKERRQ(ierr);
412
413 map<int,Gel::BlockMaterialData> &material_data = gel.blockMaterialData;
414 for(
415 map<int,Gel::BlockMaterialData>::iterator mit = material_data.begin();
416 mit!=material_data.end();
417 mit++
418 ) {
419 ierr = m_field.add_ents_to_finite_element_by_type(mit->second.tEts,MBTET,"GEL_FE"); CHKERRQ(ierr);
420 }
421
422 // Add Neumann forces, i.e. on triangles, edges and nodes.
423 ierr = MetaNeummanForces::addNeumannBCElements(m_field,"SPATIAL_POSITION"); CHKERRQ(ierr);
424 ierr = MetaNodalForces::addElement(m_field,"SPATIAL_POSITION"); CHKERRQ(ierr);
425 ierr = MetaEdgeForces::addElement(m_field,"SPATIAL_POSITION"); CHKERRQ(ierr);
426
427 // Add solvent flux element
428 {
429
430 ierr = m_field.add_finite_element("CHEMICAL_LOAD_FLUX_FE",MF_ZERO); CHKERRQ(ierr);
431 ierr = m_field.modify_finite_element_add_field_row("CHEMICAL_LOAD_FLUX_FE","CHEMICAL_LOAD"); CHKERRQ(ierr);
432 ierr = m_field.modify_finite_element_add_field_col("CHEMICAL_LOAD_FLUX_FE","CHEMICAL_LOAD"); CHKERRQ(ierr);
433 ierr = m_field.modify_finite_element_add_field_data("CHEMICAL_LOAD_FLUX_FE","CHEMICAL_LOAD"); CHKERRQ(ierr);
434 ierr = m_field.modify_finite_element_add_field_data("CHEMICAL_LOAD_FLUX_FE","MESH_NODE_POSITIONS"); CHKERRQ(ierr);
435
436 // Assume that boundary conditions are set in block containing surface
437 // triangle elements and block name is "CHEMICAL_LOAD_FLUX"
439 if(it->getName().compare(0,18,"FLUX_CHEMICAL_LOAD") == 0) {
440 vector<double> data;
441 ierr = it->getAttributes(data); CHKERRQ(ierr);
442 if(data.size()!=1) {
443 SETERRQ(PETSC_COMM_SELF,1,"Data inconsistency");
444 }
445 // Here it set how block of for heat flux is set. This is because
446 // implementation from thermal element is used to enforce this
447 // boundary condition.
448 strcpy(set_of_solvent_fluxes[it->getMeshsetId()].dAta.data.name,"HeatFlu");
449 set_of_solvent_fluxes[it->getMeshsetId()].dAta.data.flag1 = 1;
450 set_of_solvent_fluxes[it->getMeshsetId()].dAta.data.value1 = data[0];
451 //cerr << set_of_solvent_fluxes[it->getMeshsetId()].dAta << endl;
452 rval = m_field.get_moab().get_entities_by_type(
453 it->meshset,MBTRI,set_of_solvent_fluxes[it->getMeshsetId()].tRis,true
454 ); CHKERRQ_MOAB(rval);
456 set_of_solvent_fluxes[it->getMeshsetId()].tRis,MBTRI,"CHEMICAL_LOAD_FLUX_FE"
457 ); CHKERRQ(ierr);
458 }
459 }
460 }
461
462 //build finite elements
463 ierr = m_field.build_finite_elements(); CHKERRQ(ierr);
464 //build adjacencies
465 ierr = m_field.build_adjacencies(bit_level0); CHKERRQ(ierr);
466 }
467
468 }
469
470 // attach tags for each recorder
471 vector<int> tags;
472 tags.push_back(Gel::STRESSTOTAL); // ADOL-C tag used to calculate total stress
473 tags.push_back(Gel::SOLVENTFLUX);
474 tags.push_back(Gel::SOLVENTRATE);
475 tags.push_back(Gel::RESIDUALSTRAINHAT);
476
477 // Create gel instance and set operators.
478 {
479
480 gel.constitutiveEquationPtr = boost::shared_ptr<UserGelConstitutiveEquation<adouble> >(
482 );
483
484 // Set name of fields which has been choose to approximate spatial
485 // displacements, solvent concentration and internal state variables.
486 Gel::CommonData &common_data = gel.commonData;
487 common_data.spatialPositionName = "SPATIAL_POSITION";
488 common_data.spatialPositionNameDot = "SPATIAL_POSITION_DOT";
489 common_data.muName = "CHEMICAL_LOAD";
490 common_data.muNameDot = "CHEMICAL_LOAD_DOT";
491 common_data.strainHatName = "HAT_EPS";
492 common_data.strainHatNameDot = "HAT_EPS_DOT";
493
494 // Set operators to calculate field values at integration points, both for
495 // left and right hand side elements.
496 Gel::GelFE *fe_ptr[] = { &gel.feRhs, &gel.feLhs };
497 for(int ss = 0;ss<2;ss++) {
498 ierr = addHOOpsVol("MESH_NODE_POSITIONS", *fe_ptr[ss], true, true, false,
499 false);
500 CHKERRQ(ierr);
501 fe_ptr[ss]->getOpPtrVector().push_back(
502 new Gel::OpGetDataAtGaussPts("SPATIAL_POSITION",common_data,false,true)
503 );
504 fe_ptr[ss]->getOpPtrVector().push_back(
505 new Gel::OpGetDataAtGaussPts("SPATIAL_POSITION_DOT",common_data,false,true)
506 );
507 fe_ptr[ss]->getOpPtrVector().push_back(
508 new Gel::OpGetDataAtGaussPts("CHEMICAL_LOAD",common_data,true,true)
509 );
510 fe_ptr[ss]->getOpPtrVector().push_back(
511 new Gel::OpGetDataAtGaussPts("HAT_EPS",common_data,true,false,MBTET)
512 );
513 fe_ptr[ss]->getOpPtrVector().push_back(
514 new Gel::OpGetDataAtGaussPts("HAT_EPS_DOT",common_data,true,false,MBTET)
515 );
516 }
517
518
519 // Right hand side operators
520 gel.feRhs.getOpPtrVector().push_back(
521 new Gel::OpJacobian(
522 "SPATIAL_POSITION",tags,gel.constitutiveEquationPtr,gel.commonData,true,false
523 )
524 );
525 gel.feRhs.getOpPtrVector().push_back(
527 );
528 gel.feRhs.getOpPtrVector().push_back(
530 );
531 gel.feRhs.getOpPtrVector().push_back(
533 );
534 gel.feRhs.getOpPtrVector().push_back(
536 );
537
538 // Left hand side operators
539 gel.feLhs.getOpPtrVector().push_back(
540 new Gel::OpJacobian(
541 "SPATIAL_POSITION",tags,gel.constitutiveEquationPtr,gel.commonData,false,true
542 )
543 );
544 gel.feLhs.getOpPtrVector().push_back(
546 );
547 gel.feLhs.getOpPtrVector().push_back(
549 );
550 gel.feLhs.getOpPtrVector().push_back(
552 );
553 gel.feLhs.getOpPtrVector().push_back(
555 );
556 gel.feLhs.getOpPtrVector().push_back(
558 );
559 gel.feLhs.getOpPtrVector().push_back(
561 );
562 gel.feLhs.getOpPtrVector().push_back(
564 );
565
566 }
567
568 // Create discrete manager instance
569 DM dm;
570 DMType dm_name = "DMGEL";
571 {
572 ierr = DMRegister_MoFEM(dm_name); CHKERRQ(ierr);
573 ierr = DMCreate(PETSC_COMM_WORLD,&dm);CHKERRQ(ierr);
574 ierr = DMSetType(dm,dm_name);CHKERRQ(ierr);
575 ierr = DMMoFEMCreateMoFEM(dm,&m_field,dm_name,bit_level0); CHKERRQ(ierr);
576 ierr = DMSetFromOptions(dm); CHKERRQ(ierr);
577 //add elements to dm
578 ierr = DMMoFEMAddElement(dm,"GEL_FE"); CHKERRQ(ierr);
579 ierr = DMMoFEMAddElement(dm,"FORCE_FE"); CHKERRQ(ierr);
580 ierr = DMMoFEMAddElement(dm,"PRESSURE_FE"); CHKERRQ(ierr);
581 ierr = DMMoFEMAddElement(dm,"CHEMICAL_LOAD_FLUX_FE"); CHKERRQ(ierr);
582 ierr = DMSetUp(dm); CHKERRQ(ierr);
583 }
584
585 // Create matrices and vectors used for analysis
586 Vec T,F;
587 Mat A;
588 {
589 ierr = DMCreateGlobalVector_MoFEM(dm,&T); CHKERRQ(ierr);
590 ierr = VecDuplicate(T,&F); CHKERRQ(ierr);
591 ierr = DMCreateMatrix_MoFEM(dm,&A); CHKERRQ(ierr);
592 }
593
594 // Setting finite element methods for Dirichelt boundary conditions
596 m_field,"SPATIAL_POSITION",A,T,F
597 );
598 spatial_position_bc.methodsOp.push_back(new TimeForceScale("-my_displacements_history",false));
600 m_field,"CHEMICAL_LOAD","CHEMICAL_LOAD",A,T,F
601 );
602 concentration_bc.methodsOp.push_back(new TimeForceScale("-my_chemical_load_history",false));
603
604 // Setting finite element method for applying tractions
605 boost::ptr_map<string,NeummanForcesSurface> neumann_forces;
606 boost::ptr_map<string,NodalForce> nodal_forces;
607 boost::ptr_map<string,EdgeForce> edge_forces;
608 {
609 //forces and pressures on surface
610 ierr = MetaNeummanForces::setMomentumFluxOperators(m_field,neumann_forces,PETSC_NULL,"SPATIAL_POSITION"); CHKERRQ(ierr);
611 //noadl forces
612 ierr = MetaNodalForces::setOperators(m_field,nodal_forces,PETSC_NULL,"SPATIAL_POSITION"); CHKERRQ(ierr);
613 //edge forces
614 ierr = MetaEdgeForces::setOperators(m_field,edge_forces,PETSC_NULL,"SPATIAL_POSITION"); CHKERRQ(ierr);
615 for(
616 boost::ptr_map<string,NeummanForcesSurface>::iterator mit = neumann_forces.begin();
617 mit!=neumann_forces.end();mit++
618 ) {
619 mit->second->methodsOp.push_back(new TimeForceScale("-my_load_history",false));
620 }
621 for(
622 boost::ptr_map<string,NodalForce>::iterator mit = nodal_forces.begin();
623 mit!=nodal_forces.end();mit++
624 ) {
625 mit->second->methodsOp.push_back(new TimeForceScale("-my_load_history",false));
626 }
627 for(
628 boost::ptr_map<string,EdgeForce>::iterator mit = edge_forces.begin();
629 mit!=edge_forces.end();mit++
630 ) {
631 mit->second->methodsOp.push_back(new TimeForceScale("-my_load_history",false));
632 }
633 }
634
635 // Solvent surface element, flux, convection radiation
636 // TODO: add radiation and convection
637 ThermalElement::MyTriFE solvent_surface_fe(m_field);
638 {
639 map<int,ThermalElement::FluxData>::iterator sit = set_of_solvent_fluxes.begin();
640 for(;sit!=set_of_solvent_fluxes.end();sit++) {
641 // add flux operator
642 solvent_surface_fe.getOpPtrVector().push_back(
643 new ThermalElement::OpHeatFlux("CHEMICAL_LOAD",PETSC_NULL,sit->second,true)
644 );
645 }
646 }
647
648 // Add finite elements to Time Stepping Solver, using Discrete Manager Interface
649 {
650 //Rhs
651 ierr = DMMoFEMTSSetIFunction(dm,DM_NO_ELEMENT,NULL,&spatial_position_bc,NULL); CHKERRQ(ierr);
652 ierr = DMMoFEMTSSetIFunction(dm,DM_NO_ELEMENT,NULL,&concentration_bc,NULL); CHKERRQ(ierr);
653 ierr = DMMoFEMTSSetIFunction(dm,"GEL_FE",&gel.feRhs,NULL,NULL); CHKERRQ(ierr);
654 ierr = DMMoFEMTSSetIFunction(dm,"CHEMICAL_LOAD_FLUX_FE",&solvent_surface_fe,NULL,NULL); CHKERRQ(ierr);
655 {
656 boost::ptr_map<string,NeummanForcesSurface>::iterator mit = neumann_forces.begin();
657 for(;mit!=neumann_forces.end();mit++) {
658 ierr = DMMoFEMTSSetIFunction(dm,mit->first.c_str(),&mit->second->getLoopFe(),NULL,NULL); CHKERRQ(ierr);
659 }
660 }
661 {
662 boost::ptr_map<string,NodalForce>::iterator fit = nodal_forces.begin();
663 for(;fit!=nodal_forces.end();fit++) {
664 ierr = DMMoFEMTSSetIFunction(dm,fit->first.c_str(),&fit->second->getLoopFe(),NULL,NULL); CHKERRQ(ierr);
665 }
666 }
667 {
668 boost::ptr_map<string,EdgeForce>::iterator fit = edge_forces.begin();
669 for(;fit!=edge_forces.end();fit++) {
670 ierr = DMMoFEMTSSetIFunction(dm,fit->first.c_str(),&fit->second->getLoopFe(),NULL,NULL); CHKERRQ(ierr);
671 }
672 }
673 ierr = DMMoFEMTSSetIFunction(dm,DM_NO_ELEMENT,NULL,NULL,&spatial_position_bc); CHKERRQ(ierr);
674 ierr = DMMoFEMTSSetIFunction(dm,DM_NO_ELEMENT,NULL,NULL,&concentration_bc); CHKERRQ(ierr);
675
676 //Lhs
677 ierr = DMMoFEMTSSetIJacobian(dm,DM_NO_ELEMENT,NULL,&spatial_position_bc,NULL); CHKERRQ(ierr);
678 ierr = DMMoFEMTSSetIJacobian(dm,DM_NO_ELEMENT,NULL,&concentration_bc,NULL); CHKERRQ(ierr);
679 ierr = DMMoFEMTSSetIJacobian(dm,"GEL_FE",&gel.feLhs,NULL,NULL); CHKERRQ(ierr);
680 ierr = DMMoFEMTSSetIJacobian(dm,DM_NO_ELEMENT,NULL,NULL,&spatial_position_bc); CHKERRQ(ierr);
681 ierr = DMMoFEMTSSetIJacobian(dm,DM_NO_ELEMENT,NULL,NULL,&concentration_bc); CHKERRQ(ierr);
682
683 }
684
685
686 // Create Time Stepping solver
687 TS ts;
688 {
689 ierr = TSCreate(PETSC_COMM_WORLD,&ts); CHKERRQ(ierr);
690 ierr = TSSetType(ts,TSBEULER); CHKERRQ(ierr);
691 }
692
693 {
694 ierr = DMoFEMMeshToLocalVector(dm,T,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
695 ierr = DMoFEMPreProcessFiniteElements(dm,&spatial_position_bc); CHKERRQ(ierr);
696 ierr = DMoFEMPreProcessFiniteElements(dm,&concentration_bc); CHKERRQ(ierr);
697 ierr = DMoFEMMeshToGlobalVector(dm,T,INSERT_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
698 }
699
700 // Solve problem
701 {
702
703 ierr = TSSetIFunction(ts,F,PETSC_NULL,PETSC_NULL); CHKERRQ(ierr);
704 ierr = TSSetIJacobian(ts,A,A,PETSC_NULL,PETSC_NULL); CHKERRQ(ierr);
705
706 //Monitor
707 Gel::MonitorPostProc post_proc(
708 m_field,"DMGEL","GEL_FE",gel.commonData,gel.constitutiveEquationPtr,tags
709 );
710 TsCtx *ts_ctx;
712 {
713 ts_ctx->get_postProcess_to_do_Monitor().push_back(&post_proc);
714 ierr = TSMonitorSet(ts,f_TSMonitorSet,ts_ctx,PETSC_NULL); CHKERRQ(ierr);
715 }
716
717 double ftime = 1;
718 ierr = TSSetDuration(ts,PETSC_DEFAULT,ftime); CHKERRQ(ierr);
719 ierr = TSSetFromOptions(ts); CHKERRQ(ierr);
720 ierr = TSSetDM(ts,dm); CHKERRQ(ierr);
721 #if PETSC_VERSION_GE(3,7,0)
722 ierr = TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER); CHKERRQ(ierr);
723 #endif
724 ierr = TSSolve(ts,T); CHKERRQ(ierr);
725 ierr = TSGetTime(ts,&ftime); CHKERRQ(ierr);
726 PetscInt steps,snesfails,rejects,nonlinits,linits;
727 ierr = TSGetTimeStepNumber(ts,&steps); CHKERRQ(ierr);
728 ierr = TSGetSNESFailures(ts,&snesfails); CHKERRQ(ierr);
729 ierr = TSGetStepRejections(ts,&rejects); CHKERRQ(ierr);
730 ierr = TSGetSNESIterations(ts,&nonlinits); CHKERRQ(ierr);
731 ierr = TSGetKSPIterations(ts,&linits); CHKERRQ(ierr);
732 PetscPrintf(PETSC_COMM_WORLD,
733 "steps %D (%D rejected, %D SNES fails), ftime %g, nonlinits %D, linits %D\n",
734 steps,rejects,snesfails,ftime,nonlinits,linits
735 );
736 }
737
738 if(is_atom_test) {
739 double sum = 0;
740 ierr = VecSum(T,&sum); CHKERRQ(ierr);
741 ierr = PetscPrintf(PETSC_COMM_WORLD,"sum = %9.8e\n",sum); CHKERRQ(ierr);
742 double fnorm;
743 ierr = VecNorm(T,NORM_2,&fnorm); CHKERRQ(ierr);
744 ierr = PetscPrintf(PETSC_COMM_WORLD,"fnorm = %9.8e\n",fnorm); CHKERRQ(ierr);
745
746 if (fabs(sum - 6.82084809e+01) > 1e-6) {
747 SETERRQ(PETSC_COMM_WORLD,MOFEM_ATOM_TEST_INVALID,"Failed to pass test");
748 }
749 if (fabs(fnorm - 3.66184607e+01) > 1e-6) {
750 SETERRQ(PETSC_COMM_WORLD,MOFEM_ATOM_TEST_INVALID,"Failed to pass test");
751 }
752 }
753
754 // Clean and destroy
755 {
756 ierr = DMDestroy(&dm); CHKERRQ(ierr);
757 ierr = VecDestroy(&T); CHKERRQ(ierr);
758 ierr = VecDestroy(&F); CHKERRQ(ierr);
759 ierr = MatDestroy(&A); CHKERRQ(ierr);
760 //ierr = TSDestroy(&ts); CHKERRQ(ierr);
761 }
762
763 }
765
767
768 return 0;
769}
#define DM_NO_ELEMENT
Definition: DMMoFEM.hpp:10
DEPRECATED typedef DirichletSetFieldFromBlockWithFlags DirichletBCFromBlockSetFEMethodPreAndPostProc
DEPRECATED typedef DirichletSpatialPositionsBc SpatialPositionsBCFEMethodPreAndPostProc
Implementation of Gel finite element.
Post-process fields on refined mesh.
Operators and data structures for thermal analysis.
Implementation of Gel finite element.
int main()
Definition: adol-c_atom.cpp:46
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
@ MF_ZERO
Definition: definitions.h:98
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
@ L2
field with C-1 continuity
Definition: definitions.h:88
@ H1
continuous field
Definition: definitions.h:85
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:215
@ BLOCKSET
Definition: definitions.h:148
#define CHKERRQ_MOAB(a)
check error code of MoAB function
Definition: definitions.h:454
@ MOFEM_STD_EXCEPTION_THROW
Definition: definitions.h:39
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:40
constexpr int order
@ F
static char help[]
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
Definition: DMMoFEM.cpp:483
PetscErrorCode DMMoFEMTSSetIFunction(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set TS implicit function evaluation function
Definition: DMMoFEM.cpp:786
PetscErrorCode DMMoFEMCreateMoFEM(DM dm, MoFEM::Interface *m_field_ptr, const char problem_name[], const MoFEM::BitRefLevel bit_level, const MoFEM::BitRefLevel bit_mask=MoFEM::BitRefLevel().set())
Must be called by user to set MoFEM data structures.
Definition: DMMoFEM.cpp:118
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMoFEM.cpp:509
PetscErrorCode DMCreateMatrix_MoFEM(DM dm, Mat *M)
Definition: DMMoFEM.cpp:1183
PetscErrorCode DMMoFEMGetTsCtx(DM dm, MoFEM::TsCtx **ts_ctx)
get MoFEM::TsCtx data structure
Definition: DMMoFEM.cpp:1128
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:47
PetscErrorCode DMMoFEMTSSetIJacobian(DM dm, const std::string fe_name, boost::shared_ptr< MoFEM::FEMethod > method, boost::shared_ptr< MoFEM::BasicMethod > pre_only, boost::shared_ptr< MoFEM::BasicMethod > post_only)
set TS Jacobian evaluation function
Definition: DMMoFEM.cpp:839
PetscErrorCode DMCreateGlobalVector_MoFEM(DM dm, Vec *g)
DMShellSetCreateGlobalVector.
Definition: DMMoFEM.cpp:1153
PetscErrorCode DMoFEMMeshToGlobalVector(DM dm, Vec g, InsertMode mode, ScatterMode scatter_mode)
set ghosted vector values on all existing mesh entities
Definition: DMMoFEM.cpp:521
PetscErrorCode DMoFEMPreProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
Definition: DMMoFEM.cpp:532
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 bool check_field(const std::string &name) const =0
check if field is in database
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.
#define _IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet having a particular BC meshset in a moFEM field.
MoFEM::TsCtx * ts_ctx
Definition: level_set.cpp:1884
char mesh_file_name[255]
const double T
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Exceptions.hpp:74
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
MoFEMErrorCode addHOOpsVol(const std::string field, E &e, bool h1, bool hcurl, bool hdiv, bool l2)
Definition: enable_if.hpp:6
constexpr AssemblyType A
Common data for gel model.
Definition: Gels.hpp:345
definition of volume element
Definition: Gels.hpp:379
Definition: Gels.hpp:2038
Assemble matrix .
Definition: Gels.hpp:1815
Assemble matrix .
Definition: Gels.hpp:1895
Assemble matrix .
Definition: Gels.hpp:1482
Assemble matrix .
Definition: Gels.hpp:1390
Calculating right hand side.
Definition: Gels.hpp:1254
Calculate internal forces for solvent flux.
Definition: Gels.hpp:1205
Residual strain hat.
Definition: Gels.hpp:1302
Assemble internal force vector.
Definition: Gels.hpp:1154
Implementation of Gel constitutive model.
Definition: Gels.hpp:74
GelFE feLhs
Definition: Gels.hpp:436
GelFE feRhs
Definition: Gels.hpp:436
map< int, BlockMaterialData > blockMaterialData
Definition: Gels.hpp:105
@ RESIDUALSTRAINHAT
Definition: Gels.hpp:80
CommonData commonData
Definition: Gels.hpp:376
boost::shared_ptr< Gel::ConstitutiveEquation< adouble > > constitutiveEquationPtr
Definition: Gels.hpp:340
User (hackable) Gel model.
static MoFEMErrorCode setOperators(MoFEM::Interface &m_field, boost::ptr_map< std::string, EdgeForce > &edge_forces, Vec F, const std::string field_name, std::string mesh_node_positions="MESH_NODE_POSITIONS")
Set integration point operators.
Definition: EdgeForce.hpp:97
static MoFEMErrorCode addElement(MoFEM::Interface &m_field, const std::string field_name, Range *intersect_ptr=NULL)
Add element taking information from NODESET.
Definition: EdgeForce.hpp:62
static MoFEMErrorCode setOperators(MoFEM::Interface &m_field, boost::ptr_map< std::string, NodalForce > &nodal_forces, Vec F, const std::string field_name)
Set integration point operators.
Definition: NodalForce.hpp:128
static MoFEMErrorCode addElement(MoFEM::Interface &m_field, const std::string field_name, Range *intersect_ptr=NULL)
Add element taking information from NODESET.
Definition: NodalForce.hpp:92
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.
DEPRECATED MoFEMErrorCode synchronise_entities(Range &ent, int verb=DEFAULT_VERBOSITY)
DEPRECATED MoFEMErrorCode seed_ref_level_3D(const EntityHandle meshset, const BitRefLevel &bit, int verb=-1)
seed 2D entities in the meshset and their adjacencies (only TETs adjacencies) in a particular BitRefL...
boost::ptr_deque< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.
Projection of edge entities with one mid-node on hierarchical basis.
Interface for Time Stepping (TS) solver.
Definition: TsCtx.hpp:15
define surface element
operator for calculate heat flux and assemble to right hand side
Force scale operator for reading two columns.