v0.14.0
Classes | Functions | Variables
gel_analysis.cpp File Reference
User modules » Gel model

Reads cubit file and solves problem with gel material. More...

#include <MoFEM.hpp>
#include <PostProcOnRefMesh.hpp>
#include <Projection10NodeCoordsOnField.hpp>
#include <boost/numeric/ublas/vector_proxy.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/matrix_proxy.hpp>
#include <boost/numeric/ublas/vector.hpp>
#include <adolc/adolc.h>
#include <Gels.hpp>
#include <UserGelModel.hpp>
#include <boost/program_options.hpp>
#include <MethodForForceScaling.hpp>
#include <DirichletBC.hpp>
#include <SurfacePressure.hpp>
#include <NodalForce.hpp>
#include <EdgeForce.hpp>
#include <ThermalElement.hpp>
#include <TimeForceScale.hpp>
#include <boost/iostreams/tee.hpp>
#include <boost/iostreams/stream.hpp>
#include <fstream>
#include <iostream>

Go to the source code of this file.

Classes

struct  BlockData
 

Functions

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

Variables

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

Detailed Description

Reads cubit file and solves problem with gel material.

Todo:
Current version is limited only to one material. If one like to have general problem for nonlinear elasticity should implement general time dependent problem. If inertia terms need to be considered, this material should be add to nonlinear dynamics problem.
Todo:
Internal history state variables need to be statically condensed. It can be done by implementing static condensation on finite element level or by implementing pre-conditioner.

Definition in file gel_analysis.cpp.

Function Documentation

◆ main()

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

Definition at line 81 of file gel_analysis.cpp.

81  {
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",
104  mesh_file_name,
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 
249  for(_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field,BLOCKSET,it)) {
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);
319  for(_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field,BLOCKSET,it)) {
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"
438  for(_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field,BLOCKSET,it)) {
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);
455  ierr = m_field.add_ents_to_finite_element_by_type(
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> >(
481  new UserGelConstitutiveEquation<adouble>(gel.blockMaterialData)
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  fe_ptr[ss]->getOpPtrVector().push_back(
499  new Gel::OpGetDataAtGaussPts("SPATIAL_POSITION",common_data,false,true)
500  );
501  fe_ptr[ss]->getOpPtrVector().push_back(
502  new Gel::OpGetDataAtGaussPts("SPATIAL_POSITION_DOT",common_data,false,true)
503  );
504  fe_ptr[ss]->getOpPtrVector().push_back(
505  new Gel::OpGetDataAtGaussPts("CHEMICAL_LOAD",common_data,true,true)
506  );
507  fe_ptr[ss]->getOpPtrVector().push_back(
508  new Gel::OpGetDataAtGaussPts("HAT_EPS",common_data,true,false,MBTET)
509  );
510  fe_ptr[ss]->getOpPtrVector().push_back(
511  new Gel::OpGetDataAtGaussPts("HAT_EPS_DOT",common_data,true,false,MBTET)
512  );
513  }
514 
515 
516  // Right hand side operators
517  gel.feRhs.getOpPtrVector().push_back(
518  new Gel::OpJacobian(
519  "SPATIAL_POSITION",tags,gel.constitutiveEquationPtr,gel.commonData,true,false
520  )
521  );
522  gel.feRhs.getOpPtrVector().push_back(
523  new Gel::OpRhsStressTotal(gel.commonData)
524  );
525  gel.feRhs.getOpPtrVector().push_back(
526  new Gel::OpRhsSolventFlux(gel.commonData)
527  );
528  gel.feRhs.getOpPtrVector().push_back(
529  new Gel::OpRhsSolventConcetrationDot(gel.commonData)
530  );
531  gel.feRhs.getOpPtrVector().push_back(
532  new Gel::OpRhsStrainHat(gel.commonData)
533  );
534 
535  // Left hand side operators
536  gel.feLhs.getOpPtrVector().push_back(
537  new Gel::OpJacobian(
538  "SPATIAL_POSITION",tags,gel.constitutiveEquationPtr,gel.commonData,false,true
539  )
540  );
541  gel.feLhs.getOpPtrVector().push_back(
542  new Gel::OpLhsdxdx(gel.commonData)
543  );
544  gel.feLhs.getOpPtrVector().push_back(
545  new Gel::OpLhsdxdMu(gel.commonData)
546  );
547  gel.feLhs.getOpPtrVector().push_back(
548  new Gel::OpLhsdxdStrainHat(gel.commonData)
549  );
550  gel.feLhs.getOpPtrVector().push_back(
551  new Gel::OpLhsdStrainHatdStrainHat(gel.commonData)
552  );
553  gel.feLhs.getOpPtrVector().push_back(
554  new Gel::OpLhsdStrainHatdx(gel.commonData)
555  );
556  gel.feLhs.getOpPtrVector().push_back(
557  new Gel::OpLhsdMudMu(gel.commonData)
558  );
559  gel.feLhs.getOpPtrVector().push_back(
560  new Gel::OpLhsdMudx(gel.commonData)
561  );
562 
563  }
564 
565  // Create discrete manager instance
566  DM dm;
567  DMType dm_name = "DMGEL";
568  {
569  ierr = DMRegister_MoFEM(dm_name); CHKERRQ(ierr);
570  ierr = DMCreate(PETSC_COMM_WORLD,&dm);CHKERRQ(ierr);
571  ierr = DMSetType(dm,dm_name);CHKERRQ(ierr);
572  ierr = DMMoFEMCreateMoFEM(dm,&m_field,dm_name,bit_level0); CHKERRQ(ierr);
573  ierr = DMSetFromOptions(dm); CHKERRQ(ierr);
574  //add elements to dm
575  ierr = DMMoFEMAddElement(dm,"GEL_FE"); CHKERRQ(ierr);
576  ierr = DMMoFEMAddElement(dm,"FORCE_FE"); CHKERRQ(ierr);
577  ierr = DMMoFEMAddElement(dm,"PRESSURE_FE"); CHKERRQ(ierr);
578  ierr = DMMoFEMAddElement(dm,"CHEMICAL_LOAD_FLUX_FE"); CHKERRQ(ierr);
579  ierr = DMSetUp(dm); CHKERRQ(ierr);
580  }
581 
582  // Create matrices and vectors used for analysis
583  Vec T,F;
584  Mat A;
585  {
586  ierr = DMCreateGlobalVector_MoFEM(dm,&T); CHKERRQ(ierr);
587  ierr = VecDuplicate(T,&F); CHKERRQ(ierr);
588  ierr = DMCreateMatrix_MoFEM(dm,&A); CHKERRQ(ierr);
589  }
590 
591  // Setting finite element methods for Dirichelt boundary conditions
592  SpatialPositionsBCFEMethodPreAndPostProc spatial_position_bc(
593  m_field,"SPATIAL_POSITION",A,T,F
594  );
595  spatial_position_bc.methodsOp.push_back(new TimeForceScale("-my_displacements_history",false));
596  DirichletBCFromBlockSetFEMethodPreAndPostProc concentration_bc(
597  m_field,"CHEMICAL_LOAD","CHEMICAL_LOAD",A,T,F
598  );
599  concentration_bc.methodsOp.push_back(new TimeForceScale("-my_chemical_load_history",false));
600 
601  // Setting finite element method for applying tractions
602  boost::ptr_map<string,NeummanForcesSurface> neumann_forces;
603  boost::ptr_map<string,NodalForce> nodal_forces;
604  boost::ptr_map<string,EdgeForce> edge_forces;
605  {
606  //forces and pressures on surface
607  ierr = MetaNeummanForces::setMomentumFluxOperators(m_field,neumann_forces,PETSC_NULL,"SPATIAL_POSITION"); CHKERRQ(ierr);
608  //noadl forces
609  ierr = MetaNodalForces::setOperators(m_field,nodal_forces,PETSC_NULL,"SPATIAL_POSITION"); CHKERRQ(ierr);
610  //edge forces
611  ierr = MetaEdgeForces::setOperators(m_field,edge_forces,PETSC_NULL,"SPATIAL_POSITION"); CHKERRQ(ierr);
612  for(
613  boost::ptr_map<string,NeummanForcesSurface>::iterator mit = neumann_forces.begin();
614  mit!=neumann_forces.end();mit++
615  ) {
616  mit->second->methodsOp.push_back(new TimeForceScale("-my_load_history",false));
617  }
618  for(
619  boost::ptr_map<string,NodalForce>::iterator mit = nodal_forces.begin();
620  mit!=nodal_forces.end();mit++
621  ) {
622  mit->second->methodsOp.push_back(new TimeForceScale("-my_load_history",false));
623  }
624  for(
625  boost::ptr_map<string,EdgeForce>::iterator mit = edge_forces.begin();
626  mit!=edge_forces.end();mit++
627  ) {
628  mit->second->methodsOp.push_back(new TimeForceScale("-my_load_history",false));
629  }
630  }
631 
632  // Solvent surface element, flux, convection radiation
633  // TODO: add radiation and convection
634  ThermalElement::MyTriFE solvent_surface_fe(m_field);
635  {
636  map<int,ThermalElement::FluxData>::iterator sit = set_of_solvent_fluxes.begin();
637  for(;sit!=set_of_solvent_fluxes.end();sit++) {
638  // add flux operator
639  solvent_surface_fe.getOpPtrVector().push_back(
640  new ThermalElement::OpHeatFlux("CHEMICAL_LOAD",PETSC_NULL,sit->second,true)
641  );
642  }
643  }
644 
645  // Add finite elements to Time Stepping Solver, using Discrete Manager Interface
646  {
647  //Rhs
648  ierr = DMMoFEMTSSetIFunction(dm,DM_NO_ELEMENT,NULL,&spatial_position_bc,NULL); CHKERRQ(ierr);
649  ierr = DMMoFEMTSSetIFunction(dm,DM_NO_ELEMENT,NULL,&concentration_bc,NULL); CHKERRQ(ierr);
650  ierr = DMMoFEMTSSetIFunction(dm,"GEL_FE",&gel.feRhs,NULL,NULL); CHKERRQ(ierr);
651  ierr = DMMoFEMTSSetIFunction(dm,"CHEMICAL_LOAD_FLUX_FE",&solvent_surface_fe,NULL,NULL); CHKERRQ(ierr);
652  {
653  boost::ptr_map<string,NeummanForcesSurface>::iterator mit = neumann_forces.begin();
654  for(;mit!=neumann_forces.end();mit++) {
655  ierr = DMMoFEMTSSetIFunction(dm,mit->first.c_str(),&mit->second->getLoopFe(),NULL,NULL); CHKERRQ(ierr);
656  }
657  }
658  {
659  boost::ptr_map<string,NodalForce>::iterator fit = nodal_forces.begin();
660  for(;fit!=nodal_forces.end();fit++) {
661  ierr = DMMoFEMTSSetIFunction(dm,fit->first.c_str(),&fit->second->getLoopFe(),NULL,NULL); CHKERRQ(ierr);
662  }
663  }
664  {
665  boost::ptr_map<string,EdgeForce>::iterator fit = edge_forces.begin();
666  for(;fit!=edge_forces.end();fit++) {
667  ierr = DMMoFEMTSSetIFunction(dm,fit->first.c_str(),&fit->second->getLoopFe(),NULL,NULL); CHKERRQ(ierr);
668  }
669  }
670  ierr = DMMoFEMTSSetIFunction(dm,DM_NO_ELEMENT,NULL,NULL,&spatial_position_bc); CHKERRQ(ierr);
671  ierr = DMMoFEMTSSetIFunction(dm,DM_NO_ELEMENT,NULL,NULL,&concentration_bc); CHKERRQ(ierr);
672 
673  //Lhs
674  ierr = DMMoFEMTSSetIJacobian(dm,DM_NO_ELEMENT,NULL,&spatial_position_bc,NULL); CHKERRQ(ierr);
675  ierr = DMMoFEMTSSetIJacobian(dm,DM_NO_ELEMENT,NULL,&concentration_bc,NULL); CHKERRQ(ierr);
676  ierr = DMMoFEMTSSetIJacobian(dm,"GEL_FE",&gel.feLhs,NULL,NULL); CHKERRQ(ierr);
677  ierr = DMMoFEMTSSetIJacobian(dm,DM_NO_ELEMENT,NULL,NULL,&spatial_position_bc); CHKERRQ(ierr);
678  ierr = DMMoFEMTSSetIJacobian(dm,DM_NO_ELEMENT,NULL,NULL,&concentration_bc); CHKERRQ(ierr);
679 
680  }
681 
682 
683  // Create Time Stepping solver
684  TS ts;
685  {
686  ierr = TSCreate(PETSC_COMM_WORLD,&ts); CHKERRQ(ierr);
687  ierr = TSSetType(ts,TSBEULER); CHKERRQ(ierr);
688  }
689 
690  {
691  ierr = DMoFEMMeshToLocalVector(dm,T,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
692  ierr = DMoFEMPreProcessFiniteElements(dm,&spatial_position_bc); CHKERRQ(ierr);
693  ierr = DMoFEMPreProcessFiniteElements(dm,&concentration_bc); CHKERRQ(ierr);
694  ierr = DMoFEMMeshToGlobalVector(dm,T,INSERT_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
695  }
696 
697  // Solve problem
698  {
699 
700  ierr = TSSetIFunction(ts,F,PETSC_NULL,PETSC_NULL); CHKERRQ(ierr);
701  ierr = TSSetIJacobian(ts,A,A,PETSC_NULL,PETSC_NULL); CHKERRQ(ierr);
702 
703  //Monitor
704  Gel::MonitorPostProc post_proc(
705  m_field,"DMGEL","GEL_FE",gel.commonData,gel.constitutiveEquationPtr,tags
706  );
707  TsCtx *ts_ctx;
708  DMMoFEMGetTsCtx(dm,&ts_ctx);
709  {
710  ts_ctx->get_postProcess_to_do_Monitor().push_back(&post_proc);
711  ierr = TSMonitorSet(ts,f_TSMonitorSet,ts_ctx,PETSC_NULL); CHKERRQ(ierr);
712  }
713 
714  double ftime = 1;
715  ierr = TSSetDuration(ts,PETSC_DEFAULT,ftime); CHKERRQ(ierr);
716  ierr = TSSetFromOptions(ts); CHKERRQ(ierr);
717  ierr = TSSetDM(ts,dm); CHKERRQ(ierr);
718  #if PETSC_VERSION_GE(3,7,0)
719  ierr = TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER); CHKERRQ(ierr);
720  #endif
721  ierr = TSSolve(ts,T); CHKERRQ(ierr);
722  ierr = TSGetTime(ts,&ftime); CHKERRQ(ierr);
723  PetscInt steps,snesfails,rejects,nonlinits,linits;
724  ierr = TSGetTimeStepNumber(ts,&steps); CHKERRQ(ierr);
725  ierr = TSGetSNESFailures(ts,&snesfails); CHKERRQ(ierr);
726  ierr = TSGetStepRejections(ts,&rejects); CHKERRQ(ierr);
727  ierr = TSGetSNESIterations(ts,&nonlinits); CHKERRQ(ierr);
728  ierr = TSGetKSPIterations(ts,&linits); CHKERRQ(ierr);
729  PetscPrintf(PETSC_COMM_WORLD,
730  "steps %D (%D rejected, %D SNES fails), ftime %g, nonlinits %D, linits %D\n",
731  steps,rejects,snesfails,ftime,nonlinits,linits
732  );
733  }
734 
735  if(is_atom_test) {
736  double sum = 0;
737  ierr = VecSum(T,&sum); CHKERRQ(ierr);
738  ierr = PetscPrintf(PETSC_COMM_WORLD,"sum = %9.8e\n",sum); CHKERRQ(ierr);
739  double fnorm;
740  ierr = VecNorm(T,NORM_2,&fnorm); CHKERRQ(ierr);
741  ierr = PetscPrintf(PETSC_COMM_WORLD,"fnorm = %9.8e\n",fnorm); CHKERRQ(ierr);
742 
743  if(fabs(sum-7.57417437e+01)>1e-6) {
744  SETERRQ(PETSC_COMM_WORLD,MOFEM_ATOM_TEST_INVALID,"Failed to pass test");
745  }
746  if(fabs(fnorm-3.67577050e+01)>1e-6) {
747  SETERRQ(PETSC_COMM_WORLD,MOFEM_ATOM_TEST_INVALID,"Failed to pass test");
748  }
749  }
750 
751  // Clean and destroy
752  {
753  ierr = DMDestroy(&dm); CHKERRQ(ierr);
754  ierr = VecDestroy(&T); CHKERRQ(ierr);
755  ierr = VecDestroy(&F); CHKERRQ(ierr);
756  ierr = MatDestroy(&A); CHKERRQ(ierr);
757  //ierr = TSDestroy(&ts); CHKERRQ(ierr);
758  }
759 
760  }
761  CATCH_ERRORS;
762 
763  MoFEM::Core::Finalize();
764 
765  return 0;
766 }

Variable Documentation

◆ help

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

Definition at line 72 of file gel_analysis.cpp.

GelModule::UserGelConstitutiveEquation
User (hackable) Gel model.
Definition: UserGelModel.hpp:39
MYPCOMM_INDEX
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:228
MoFEM::CoreInterface::loop_dofs
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.
GelModule::Gel::OpLhsdStrainHatdx
Assemble matrix .
Definition: Gels.hpp:1753
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:82
H1
@ H1
continuous field
Definition: definitions.h:85
MOFEM_STD_EXCEPTION_THROW
@ MOFEM_STD_EXCEPTION_THROW
Definition: definitions.h:39
GelModule::Gel::OpJacobian
Definition: Gels.hpp:541
EntityHandle
MetaNodalForces::addElement
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
CHKERRQ_MOAB
#define CHKERRQ_MOAB(a)
check error code of MoAB function
Definition: definitions.h:467
GelModule::Gel::CommonData::spatialPositionName
string spatialPositionName
Definition: Gels.hpp:347
MoFEM::CoreInterface::modify_finite_element_add_field_row
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
L2
@ L2
field with C-1 continuity
Definition: definitions.h:88
DirichletBCFromBlockSetFEMethodPreAndPostProc
DEPRECATED typedef DirichletSetFieldFromBlockWithFlags DirichletBCFromBlockSetFEMethodPreAndPostProc
Definition: DirichletBC.hpp:381
ThermalElement::OpHeatFlux
operator for calculate heat flux and assemble to right hand side
Definition: ThermalElement.hpp:380
A
constexpr AssemblyType A
Definition: operators_tests.cpp:30
MoFEM::DMoFEMMeshToLocalVector
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:523
MoFEM::Projection10NodeCoordsOnField
Projection of edge entities with one mid-node on hierarchical basis.
Definition: Projection10NodeCoordsOnField.hpp:24
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
GelModule::Gel::OpRhsStrainHat
Residual strain hat.
Definition: Gels.hpp:1311
ts_ctx
MoFEM::TsCtx * ts_ctx
Definition: level_set.cpp:1932
GelModule::Gel::OpLhsdxdMu
Assemble matrix .
Definition: Gels.hpp:1497
MoFEM::CoreInterface::add_ents_to_field_by_type
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.
MoFEM::DMMoFEMAddElement
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
Definition: DMMoFEM.cpp:497
order
constexpr int order
Definition: dg_projection.cpp:18
GelModule::Gel
Implementation of Gel constitutive model.
Definition: Gels.hpp:74
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:2002
MetaNodalForces::setOperators
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
GelModule::Gel::OpLhsdMudx
Assemble matrix .
Definition: Gels.hpp:1925
MoFEM::DMCreateGlobalVector_MoFEM
PetscErrorCode DMCreateGlobalVector_MoFEM(DM dm, Vec *g)
DMShellSetCreateGlobalVector.
Definition: DMMoFEM.cpp:1167
MoFEM::CoreInterface::add_ents_to_finite_element_by_type
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
MoFEM::Exceptions::rval
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Exceptions.hpp:74
MoFEM::DMCreateMatrix_MoFEM
PetscErrorCode DMCreateMatrix_MoFEM(DM dm, Mat *M)
Definition: DMMoFEM.cpp:1197
MoFEM::CoreInterface::add_finite_element
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
MoFEM::DMMoFEMTSSetIJacobian
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:853
MoFEM::CoreInterface::modify_finite_element_add_field_col
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
TimeForceScale
Force scale operator for reading two columns.
Definition: TimeForceScale.hpp:18
MoFEM::CoreInterface::build_finite_elements
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
MoFEM::CoreInterface::add_field
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.
GelModule::Gel::MonitorPostProc
Definition: Gels.hpp:2071
GelModule::Gel::CommonData
Common data for gel model.
Definition: Gels.hpp:345
DM_NO_ELEMENT
#define DM_NO_ELEMENT
Definition: DMMoFEM.hpp:10
MoFEM::DMoFEMMeshToGlobalVector
PetscErrorCode DMoFEMMeshToGlobalVector(DM dm, Vec g, InsertMode mode, ScatterMode scatter_mode)
set ghosted vector values on all existing mesh entities
Definition: DMMoFEM.cpp:535
MoFEM::TsCtx
Interface for Time Stepping (TS) solver.
Definition: TsCtx.hpp:17
MoFEM::DMRegister_MoFEM
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:43
ThermalElement::MyTriFE
define surface element
Definition: ThermalElement.hpp:60
MoFEM::DeprecatedCoreInterface::synchronise_entities
DEPRECATED MoFEMErrorCode synchronise_entities(Range &ent, int verb=DEFAULT_VERBOSITY)
Definition: DeprecatedCoreInterface.cpp:41
MoFEM::DMoFEMPreProcessFiniteElements
PetscErrorCode DMoFEMPreProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
Definition: DMMoFEM.cpp:546
GelModule::Gel::OpGetDataAtGaussPts
Definition: Gels.hpp:444
MoFEM::CoreInterface::check_field
virtual bool check_field(const std::string &name) const =0
check if field is in database
MoFEM::DMMoFEMCreateMoFEM
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:114
GelModule::Gel::OpLhsdMudMu
Assemble matrix .
Definition: Gels.hpp:1842
mesh_file_name
char mesh_file_name[255]
Definition: mesh_smoothing.cpp:23
MetaEdgeForces::addElement
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
Range
GelModule::Gel::OpRhsSolventConcetrationDot
Calculating right hand side.
Definition: Gels.hpp:1260
MoFEM::CoreTmp< 0 >::Initialize
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
MF_ZERO
@ MF_ZERO
Definition: definitions.h:111
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:385
_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_
#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.
Definition: MeshsetsManager.hpp:71
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1140
BLOCKSET
@ BLOCKSET
Definition: definitions.h:161
GelModule::Gel::OpLhsdxdStrainHat
Assemble matrix .
Definition: Gels.hpp:1580
MoFEM::Exceptions::ierr
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
SpatialPositionsBCFEMethodPreAndPostProc
DEPRECATED typedef DirichletSpatialPositionsBc SpatialPositionsBCFEMethodPreAndPostProc
Definition: DirichletBC.hpp:237
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
MoFEM::DMMoFEMTSSetIFunction
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:800
MoFEM::ForcesAndSourcesCore::getOpPtrVector
boost::ptr_deque< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.
Definition: ForcesAndSourcesCore.hpp:83
EigenMatrix::Vec
const FTensor::Tensor2< T, Dim, Dim > Vec
Definition: MatrixFunction.hpp:66
MoFEM::CoreInterface::build_fields
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
MoFEM::DeprecatedCoreInterface::seed_ref_level_3D
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...
Definition: DeprecatedCoreInterface.cpp:18
GelModule::Gel::OpRhsStressTotal
Assemble internal force vector.
Definition: Gels.hpp:1154
MoFEM::CoreInterface::modify_finite_element_add_field_data
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string name_filed)=0
set finite element field data
MoFEM::DMMoFEMGetTsCtx
PetscErrorCode DMMoFEMGetTsCtx(DM dm, MoFEM::TsCtx **ts_ctx)
get MoFEM::TsCtx data structure
Definition: DMMoFEM.cpp:1142
GelModule::Gel::OpLhsdStrainHatdStrainHat
Assemble matrix .
Definition: Gels.hpp:1671
MoFEM::Types::BitRefLevel
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:40
MOFEM_ATOM_TEST_INVALID
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:40
MoFEM::CoreInterface::build_adjacencies
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
help
static char help[]
Definition: gel_analysis.cpp:72
MoFEM::CoreInterface::set_field_order
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.
GelModule::Gel::OpRhsSolventFlux
Calculate internal forces for solvent flux.
Definition: Gels.hpp:1208
GelModule::Gel::OpLhsdxdx
Assemble matrix .
Definition: Gels.hpp:1402
GelModule::Gel::GelFE
definition of volume element
Definition: Gels.hpp:379
F
@ F
Definition: free_surface.cpp:394
MetaEdgeForces::setOperators
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
Generated by Doxygen 1.8.17 and hosted at University of Glasgow