v0.6.9
Classes | Functions | Variables
gel_analysis.cpp File Reference

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.

Minimal surface area

Implementation is based on: https://www.dealii.org/developer/doxygen/deal.II/step_15.html Look for formulation details on deal.II web page.

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.
Todo:

Make this example with heterogenous approx. apace

Make this example with several refined meshes

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  PetscInitialize(&argc,&argv,(char *)0,help);
84 
85  moab::Core mb_instance;
86  moab::Interface& moab = mb_instance;
87 
88  PetscBool flg_gel_config,flg_file;
89  char mesh_file_name[255];
90  char gel_config_file[255];
91  int order = 2;
92  PetscBool is_partitioned = PETSC_FALSE;
93  PetscBool is_atom_test = PETSC_FALSE;
94 
95  {
96 
97  ierr = PetscOptionsBegin(PETSC_COMM_WORLD,"","Elastic Config","none"); CHKERRQ(ierr);
98  ierr = PetscOptionsString(
99  "-my_file",
100  "mesh file name","",
101  "mesh.h5m",
102  mesh_file_name,
103  255,
104  &flg_file
105  ); CHKERRQ(ierr);
106  ierr = PetscOptionsInt(
107  "-my_order",
108  "default approximation order",
109  "",
110  order,
111  &order,
112  PETSC_NULL
113  ); CHKERRQ(ierr);
114 
115  ierr = PetscOptionsBool(
116  "-my_is_partitioned",
117  "set if mesh is partitioned (this result that each process keeps only part of the mes","",
118  PETSC_FALSE,&is_partitioned,PETSC_NULL
119  ); CHKERRQ(ierr);
120  ierr = PetscOptionsString(
121  "-my_gel_config",
122  "gel configuration file name","",
123  "gel_config.in",gel_config_file,255,&flg_gel_config
124  ); CHKERRQ(ierr);
125 
126  ierr = PetscOptionsBool(
127  "-my_is_atom_test",
128  "is used with testing, exit with error when diverged","",
129  PETSC_FALSE,&is_atom_test,PETSC_NULL
130  ); CHKERRQ(ierr);
131 
132  ierr = PetscOptionsEnd(); CHKERRQ(ierr);
133 
134  //Reade parameters from line command
135  if(flg_file != PETSC_TRUE) {
136  SETERRQ(PETSC_COMM_SELF,1,"*** ERROR -my_file (MESH FILE NEEDED)");
137  }
138 
139  if(is_partitioned == PETSC_TRUE) {
140  //Read mesh to MOAB
141  const char *option;
142  option =
143  "PARALLEL=BCAST_DELETE;"
144  "PARALLEL_RESOLVE_SHARED_ENTS;"
145  "PARTITION=PARALLEL_PARTITION;";
146  rval = moab.load_file(mesh_file_name, 0, option); CHKERRQ_MOAB(rval);
147  } else {
148  const char *option;
149  option = "";
150  rval = moab.load_file(mesh_file_name, 0, option); CHKERRQ_MOAB(rval);
151  }
152  ParallelComm* pcomm = ParallelComm::get_pcomm(&moab,MYPCOMM_INDEX);
153  if(pcomm == NULL) pcomm = new ParallelComm(&moab,PETSC_COMM_WORLD);
154  }
155 
156  MoFEM::Core core(moab);
157  MoFEM::Interface& m_field = core;
158 
159  // Seed all mesh entities to MoFEM database, those entities can be potentially used as finite elements
160  // or as entities which carry some approximation field.
161  BitRefLevel bit_level0;
162  bit_level0.set(0);
163  ierr = m_field.seed_ref_level_3D(0,bit_level0); CHKERRQ(ierr);
164 
165  // Define fields and finite elements
166  Gel gel(m_field);
167  map<int,ThermalElement::FluxData> set_of_solvent_fluxes;
168  {
169 
170  // Set approximation fields
171  {
172 
173  // Add fields
174  bool check_if_spatial_field_exist = m_field.check_field("SPATIAL_POSITION");
175  ierr = m_field.add_field("SPATIAL_POSITION",H1,AINSWORTH_LEGENDRE_BASE,3,MB_TAG_SPARSE,MF_ZERO); CHKERRQ(ierr);
176  ierr = m_field.add_field("CHEMICAL_LOAD",H1,AINSWORTH_LEGENDRE_BASE,1,MB_TAG_SPARSE,MF_ZERO); CHKERRQ(ierr);
177  ierr = m_field.add_field("HAT_EPS",L2,AINSWORTH_LEGENDRE_BASE,6,MB_TAG_SPARSE,MF_ZERO); CHKERRQ(ierr);
178  ierr = m_field.add_field("SPATIAL_POSITION_DOT",H1,AINSWORTH_LEGENDRE_BASE,3,MB_TAG_SPARSE,MF_ZERO); CHKERRQ(ierr);
179  ierr = m_field.add_field("CHEMICAL_LOAD_DOT",H1,AINSWORTH_LEGENDRE_BASE,1,MB_TAG_SPARSE,MF_ZERO); CHKERRQ(ierr);
180  ierr = m_field.add_field("HAT_EPS_DOT",L2,AINSWORTH_LEGENDRE_BASE,6,MB_TAG_SPARSE,MF_ZERO); CHKERRQ(ierr);
181 
182  //Add field H1 space rank 3 to approximate geometry using hierarchical basis
183  //For 10 node tetrahedral, before use, geometry is projected on that field (see below)
184  ierr = m_field.add_field("MESH_NODE_POSITIONS",H1,AINSWORTH_LEGENDRE_BASE,3,MB_TAG_SPARSE,MF_ZERO); CHKERRQ(ierr);
185 
186  //meshset consisting all entities in mesh
187  EntityHandle root_set = moab.get_root_set();
188  //Add entities to field (root_mesh, i.e. on all mesh entities fields are approx.)
189  ierr = m_field.add_ents_to_field_by_type(root_set,MBTET,"SPATIAL_POSITION"); CHKERRQ(ierr);
190  ierr = m_field.add_ents_to_field_by_type(root_set,MBTET,"CHEMICAL_LOAD"); CHKERRQ(ierr);
191  ierr = m_field.add_ents_to_field_by_type(root_set,MBTET,"HAT_EPS"); CHKERRQ(ierr);
192  ierr = m_field.add_ents_to_field_by_type(root_set,MBTET,"SPATIAL_POSITION_DOT"); CHKERRQ(ierr);
193  ierr = m_field.add_ents_to_field_by_type(root_set,MBTET,"CHEMICAL_LOAD_DOT"); CHKERRQ(ierr);
194  ierr = m_field.add_ents_to_field_by_type(root_set,MBTET,"HAT_EPS_DOT"); CHKERRQ(ierr);
195 
196  // Set approximation order. Solvent concentration has one order less than
197  // order of of spatial position field. Tests need to be maid if that is
198  // good enough to have stability for this type of problem. If not bubble
199  // functions could be easily added by increasing approximate order for
200  // volume.
201 
202  ierr = m_field.set_field_order(root_set,MBTET,"SPATIAL_POSITION",order,2); CHKERRQ(ierr);
203  ierr = m_field.set_field_order(root_set,MBTRI,"SPATIAL_POSITION",order,2); CHKERRQ(ierr);
204  ierr = m_field.set_field_order(root_set,MBEDGE,"SPATIAL_POSITION",order,2); CHKERRQ(ierr);
205  ierr = m_field.set_field_order(root_set,MBVERTEX,"SPATIAL_POSITION",1); CHKERRQ(ierr);
206 
207  ierr = m_field.set_field_order(root_set,MBTET,"SPATIAL_POSITION_DOT",order); CHKERRQ(ierr);
208  ierr = m_field.set_field_order(root_set,MBTRI,"SPATIAL_POSITION_DOT",order); CHKERRQ(ierr);
209  ierr = m_field.set_field_order(root_set,MBEDGE,"SPATIAL_POSITION_DOT",order); CHKERRQ(ierr);
210  ierr = m_field.set_field_order(root_set,MBVERTEX,"SPATIAL_POSITION_DOT",1); CHKERRQ(ierr);
211 
212  ierr = m_field.set_field_order(root_set,MBTET,"CHEMICAL_LOAD",order-1); CHKERRQ(ierr);
213  ierr = m_field.set_field_order(root_set,MBTRI,"CHEMICAL_LOAD",order-1); CHKERRQ(ierr);
214  ierr = m_field.set_field_order(root_set,MBEDGE,"CHEMICAL_LOAD",order-1); CHKERRQ(ierr);
215  ierr = m_field.set_field_order(root_set,MBVERTEX,"CHEMICAL_LOAD",1); CHKERRQ(ierr);
216 
217  ierr = m_field.set_field_order(root_set,MBTET,"CHEMICAL_LOAD_DOT",order-1); CHKERRQ(ierr);
218  ierr = m_field.set_field_order(root_set,MBTRI,"CHEMICAL_LOAD_DOT",order-1); CHKERRQ(ierr);
219  ierr = m_field.set_field_order(root_set,MBEDGE,"CHEMICAL_LOAD_DOT",order-1); CHKERRQ(ierr);
220  ierr = m_field.set_field_order(root_set,MBVERTEX,"CHEMICAL_LOAD_DOT",1); CHKERRQ(ierr);
221 
222  ierr = m_field.set_field_order(root_set,MBTET,"HAT_EPS",order-1); CHKERRQ(ierr);
223  ierr = m_field.set_field_order(root_set,MBTET,"HAT_EPS_DOT",order-1); CHKERRQ(ierr);
224 
225  //gemetry approximation is set to 2nd oreder
226  ierr = m_field.add_ents_to_field_by_type(root_set,MBTET,"MESH_NODE_POSITIONS"); CHKERRQ(ierr);
227  ierr = m_field.set_field_order(0,MBTET,"MESH_NODE_POSITIONS",2); CHKERRQ(ierr);
228  ierr = m_field.set_field_order(0,MBTRI,"MESH_NODE_POSITIONS",2); CHKERRQ(ierr);
229  ierr = m_field.set_field_order(0,MBEDGE,"MESH_NODE_POSITIONS",2); CHKERRQ(ierr);
230  ierr = m_field.set_field_order(0,MBVERTEX,"MESH_NODE_POSITIONS",1); CHKERRQ(ierr);
231 
232  try {
233 
234  map<int,BlockData> block_data;
235  map<int,Gel::BlockMaterialData> &material_data = gel.blockMaterialData;
236 
237  ifstream ini_file(gel_config_file);
238  po::variables_map vm;
239  po::options_description config_file_options;
240 
241  int def_block = 0;
242  Tag th_block_id;
243  rval = moab.tag_get_handle(
244  "BLOCK_ID",1,MB_TYPE_INTEGER,th_block_id,MB_TAG_CREAT|MB_TAG_SPARSE,&def_block
245  ); CHKERRQ_MOAB(rval);
246 
248 
249  const string &name = it->getName();
250  if(name.compare(0,3,"GEL")!=0) continue;
251 
252  int block_id = it->getMeshsetId();
253 
254  ostringstream str_order;
255  str_order << "block_" << block_id << ".oRder";
256  config_file_options.add_options()
257  (str_order.str().c_str(),po::value<int>(
258  &block_data[block_id].oRder
259  )->default_value(order));
260  Range block_ents;
261  rval = moab.get_entities_by_handle(
262  it->getMeshset(),block_ents,true
263  ); CHKERRQ_MOAB(rval);
264  material_data[block_id].tEts = block_ents.subset_by_type(MBTET);
265  std::vector<int> block_id_vec(material_data[block_id].tEts.size(),it->getMeshsetId());
266  // cerr << "AAA " << block_id_vec.size() << " " << material_data[block_id].tEts.size() << endl;
267  rval = moab.tag_set_data(
268  th_block_id,material_data[block_id].tEts,&*block_id_vec.begin()
269  );
270 
271  ostringstream str_g_alpha;
272  str_g_alpha << "block_" << block_id << ".gAlpha";
273  config_file_options.add_options()
274  (str_g_alpha.str().c_str(),po::value<double>(&material_data[block_id].gAlpha)->default_value(1));
275  ostringstream str_v_alpha;
276  str_v_alpha << "block_" << block_id << ".vAlpha";
277  config_file_options.add_options()
278  (str_v_alpha.str().c_str(),po::value<double>(&material_data[block_id].vAlpha)->default_value(0));
279  ostringstream str_g_beta;
280  str_g_beta << "block_" << block_id << ".gBeta";
281  config_file_options.add_options()
282  (str_g_beta.str().c_str(),po::value<double>(&material_data[block_id].gBeta)->default_value(1));
283  ostringstream str_v_beta;
284  str_v_beta << "block_" << block_id << ".vBeta";
285  config_file_options.add_options()
286  (str_v_beta.str().c_str(),po::value<double>(&material_data[block_id].vBeta)->default_value(0));
287  ostringstream str_g_beta_hat;
288  str_g_beta_hat << "block_" << block_id << ".gBetaHat";
289  config_file_options.add_options()
290  (str_g_beta_hat.str().c_str(),po::value<double>(&material_data[block_id].gBetaHat)->default_value(1));
291  ostringstream str_v_beta_hat;
292  str_v_beta_hat << "block_" << block_id << ".vBetaHat";
293  config_file_options.add_options()
294  (str_v_beta_hat.str().c_str(),po::value<double>(&material_data[block_id].vBetaHat)->default_value(0));
295  ostringstream str_omega;
296  str_omega << "block_" << block_id << ".oMega";
297  config_file_options.add_options()
298  (str_omega.str().c_str(),po::value<double>(&material_data[block_id].oMega)->default_value(1));
299  ostringstream str_viscosity;
300  str_viscosity << "block_" << block_id << ".vIscosity";
301  config_file_options.add_options()
302  (str_viscosity.str().c_str(),po::value<double>(&material_data[block_id].vIscosity)->default_value(1));
303  ostringstream str_permeability;
304  str_permeability << "block_" << block_id << ".pErmeability";
305  config_file_options.add_options()
306  (str_permeability.str().c_str(),po::value<double>(&material_data[block_id].pErmeability)->default_value(1));
307  ostringstream str_mu0;
308  str_mu0 << "block_" << block_id << ".mU0";
309  config_file_options.add_options()
310  (str_mu0.str().c_str(),po::value<double>(&material_data[block_id].mU0)->default_value(1));
311 
312  }
313 
314  po::parsed_options parsed = parse_config_file(ini_file,config_file_options,true);
315  store(parsed,vm);
316  po::notify(vm);
318 
319  const string &name = it->getName();
320  if(name.compare(0,3,"GEL")!=0) continue;
321 
322  if(block_data[it->getMeshsetId()].oRder == -1) continue;
323  if(block_data[it->getMeshsetId()].oRder == order) continue;
324  PetscPrintf(PETSC_COMM_WORLD,"Set block %d order to %d\n",it->getMeshsetId(),block_data[it->getMeshsetId()].oRder);
325  Range block_ents;
326  rval = moab.get_entities_by_handle(it->getMeshset(),block_ents,true); CHKERRQ_MOAB(rval);
327  Range ents_to_set_order;
328  rval = moab.get_adjacencies(block_ents,3,false,ents_to_set_order,moab::Interface::UNION); CHKERRQ_MOAB(rval);
329  ents_to_set_order = ents_to_set_order.subset_by_type(MBTET);
330  rval = moab.get_adjacencies(block_ents,2,false,ents_to_set_order,moab::Interface::UNION); CHKERRQ_MOAB(rval);
331  rval = moab.get_adjacencies(block_ents,1,false,ents_to_set_order,moab::Interface::UNION); CHKERRQ_MOAB(rval);
332  ierr = m_field.synchronise_entities(ents_to_set_order); CHKERRQ(ierr);
333  ierr = m_field.set_field_order(ents_to_set_order,"SPATIAL_POSITION",block_data[it->getMeshsetId()].oRder); CHKERRQ(ierr);
334  ierr = m_field.set_field_order(ents_to_set_order,"CHEMICAL_LOAD",block_data[it->getMeshsetId()].oRder-1); CHKERRQ(ierr);
335  ierr = m_field.set_field_order(ents_to_set_order,"HAT_EPS",block_data[it->getMeshsetId()].oRder-1); CHKERRQ(ierr);
336  ierr = m_field.set_field_order(ents_to_set_order,"SPATIAL_POSITION_DOT",block_data[it->getMeshsetId()].oRder); CHKERRQ(ierr);
337  ierr = m_field.set_field_order(ents_to_set_order,"CHEMICAL_LOAD_DOT",block_data[it->getMeshsetId()].oRder-1); CHKERRQ(ierr);
338  ierr = m_field.set_field_order(ents_to_set_order,"HAT_EPS_DOT",block_data[it->getMeshsetId()].oRder-1); CHKERRQ(ierr);
339  }
340  vector<string> additional_parameters;
341  additional_parameters = collect_unrecognized(parsed.options,po::include_positional);
342  for(vector<string>::iterator vit = additional_parameters.begin();
343  vit!=additional_parameters.end();vit++) {
344  ierr = PetscPrintf(PETSC_COMM_WORLD,"** WARNING Unrecognised option %s\n",vit->c_str()); CHKERRQ(ierr);
345  }
346  } catch (const std::exception& ex) {
347  ostringstream ss;
348  ss << ex.what() << endl;
349  SETERRQ(PETSC_COMM_SELF,MOFEM_STD_EXCEPTION_THROW,ss.str().c_str());
350  }
351 
352  for(
353  map<int,Gel::BlockMaterialData>::iterator mit = gel.blockMaterialData.begin();
354  mit!=gel.blockMaterialData.end();mit++
355  ) {
356  ierr = PetscPrintf(PETSC_COMM_WORLD,"** Block Id %d\n",mit->first); CHKERRQ(ierr);
357  PetscPrintf(PETSC_COMM_WORLD,"vAlpha = %5.4g\n",mit->second.vAlpha);
358  PetscPrintf(PETSC_COMM_WORLD,"gAlpha = %5.4g\n",mit->second.gAlpha);
359  PetscPrintf(PETSC_COMM_WORLD,"vBeta = %5.4g\n",mit->second.vBeta);
360  PetscPrintf(PETSC_COMM_WORLD,"gBeta = %5.4g\n",mit->second.gBeta);
361  PetscPrintf(PETSC_COMM_WORLD,"vBetaHat = %5.4g\n",mit->second.vBetaHat);
362  PetscPrintf(PETSC_COMM_WORLD,"gBetaHat = %5.4g\n",mit->second.gBetaHat);
363  PetscPrintf(PETSC_COMM_WORLD,"pErmeability = %5.4g\n",mit->second.pErmeability);
364  PetscPrintf(PETSC_COMM_WORLD,"vIscosity = %5.4g\n",mit->second.vIscosity);
365  PetscPrintf(PETSC_COMM_WORLD,"oMega = %5.4g\n",mit->second.oMega);
366  PetscPrintf(PETSC_COMM_WORLD,"mU0 = %5.4g\n",mit->second.mU0);
367  }
368 
369  ierr = m_field.build_fields(); CHKERRQ(ierr);
370 
371  // Sett geometry approximation and initial spatial positions
372  // 10 node tets
373  if (!check_if_spatial_field_exist) {
374  Projection10NodeCoordsOnField ent_method_material(m_field, "MESH_NODE_POSITIONS");
375  ierr = m_field.loop_dofs("MESH_NODE_POSITIONS", ent_method_material); CHKERRQ(ierr);
376  Projection10NodeCoordsOnField ent_method_spatial(m_field, "SPATIAL_POSITION");
377  ierr = m_field.loop_dofs("SPATIAL_POSITION", ent_method_spatial); CHKERRQ(ierr);
378  // Set value to chemical load
379  map<int,Gel::BlockMaterialData> &material_data = gel.blockMaterialData;
380  for(
381  map<int,Gel::BlockMaterialData>::iterator mit = material_data.begin();
382  mit!=material_data.end();mit++
383  ) {
384  // FIXME mu0 at border in undefined, if two blackset have different mu0
385  Range vertices;
386  rval = moab.get_connectivity(mit->second.tEts,vertices,true); CHKERRQ_MOAB(rval);
387  ierr = m_field.set_field(mit->second.mU0,MBVERTEX,vertices,"CHEMICAL_LOAD"); CHKERRQ(ierr);
388  }
389  }
390 
391  }
392 
393  //Set finite elements. The primary element is GEL_FE in addition elements
394  //for applying tractions and fluxes are added.
395  {
396  ierr = m_field.add_finite_element("GEL_FE",MF_ZERO); CHKERRQ(ierr);
397  ierr = m_field.modify_finite_element_add_field_row("GEL_FE","SPATIAL_POSITION"); CHKERRQ(ierr);
398  ierr = m_field.modify_finite_element_add_field_col("GEL_FE","SPATIAL_POSITION"); CHKERRQ(ierr);
399  ierr = m_field.modify_finite_element_add_field_row("GEL_FE","CHEMICAL_LOAD"); CHKERRQ(ierr);
400  ierr = m_field.modify_finite_element_add_field_col("GEL_FE","CHEMICAL_LOAD"); CHKERRQ(ierr);
401  ierr = m_field.modify_finite_element_add_field_row("GEL_FE","HAT_EPS"); CHKERRQ(ierr);
402  ierr = m_field.modify_finite_element_add_field_col("GEL_FE","HAT_EPS"); CHKERRQ(ierr);
403  ierr = m_field.modify_finite_element_add_field_data("GEL_FE","SPATIAL_POSITION"); CHKERRQ(ierr);
404  ierr = m_field.modify_finite_element_add_field_data("GEL_FE","SPATIAL_POSITION_DOT"); CHKERRQ(ierr);
405  ierr = m_field.modify_finite_element_add_field_data("GEL_FE","CHEMICAL_LOAD"); CHKERRQ(ierr);
406  ierr = m_field.modify_finite_element_add_field_data("GEL_FE","CHEMICAL_LOAD_DOT"); CHKERRQ(ierr);
407  ierr = m_field.modify_finite_element_add_field_data("GEL_FE","HAT_EPS"); CHKERRQ(ierr);
408  ierr = m_field.modify_finite_element_add_field_data("GEL_FE","HAT_EPS_DOT"); CHKERRQ(ierr);
409  ierr = m_field.modify_finite_element_add_field_data("GEL_FE","MESH_NODE_POSITIONS"); CHKERRQ(ierr);
410 
411  map<int,Gel::BlockMaterialData> &material_data = gel.blockMaterialData;
412  for(
413  map<int,Gel::BlockMaterialData>::iterator mit = material_data.begin();
414  mit!=material_data.end();
415  mit++
416  ) {
417  ierr = m_field.add_ents_to_finite_element_by_type(mit->second.tEts,MBTET,"GEL_FE"); CHKERRQ(ierr);
418  }
419 
420  // Add Neumann forces, i.e. on triangles, edges and nodes.
421  ierr = MetaNeummanForces::addNeumannBCElements(m_field,"SPATIAL_POSITION"); CHKERRQ(ierr);
422  ierr = MetaNodalForces::addElement(m_field,"SPATIAL_POSITION"); CHKERRQ(ierr);
423  ierr = MetaEdgeForces::addElement(m_field,"SPATIAL_POSITION"); CHKERRQ(ierr);
424 
425  // Add solvent flux element
426  {
427 
428  ierr = m_field.add_finite_element("CHEMICAL_LOAD_FLUX_FE",MF_ZERO); CHKERRQ(ierr);
429  ierr = m_field.modify_finite_element_add_field_row("CHEMICAL_LOAD_FLUX_FE","CHEMICAL_LOAD"); CHKERRQ(ierr);
430  ierr = m_field.modify_finite_element_add_field_col("CHEMICAL_LOAD_FLUX_FE","CHEMICAL_LOAD"); CHKERRQ(ierr);
431  ierr = m_field.modify_finite_element_add_field_data("CHEMICAL_LOAD_FLUX_FE","CHEMICAL_LOAD"); CHKERRQ(ierr);
432  ierr = m_field.modify_finite_element_add_field_data("CHEMICAL_LOAD_FLUX_FE","MESH_NODE_POSITIONS"); CHKERRQ(ierr);
433 
434  // Assume that boundary conditions are set in block containing surface
435  // triangle elements and block name is "CHEMICAL_LOAD_FLUX"
437  if(it->getName().compare(0,18,"FLUX_CHEMICAL_LOAD") == 0) {
438  vector<double> data;
439  ierr = it->getAttributes(data); CHKERRQ(ierr);
440  if(data.size()!=1) {
441  SETERRQ(PETSC_COMM_SELF,1,"Data inconsistency");
442  }
443  // Here it set how block of for heat flux is set. This is because
444  // implementation from thermal element is used to enforce this
445  // boundary condition.
446  strcpy(set_of_solvent_fluxes[it->getMeshsetId()].dAta.data.name,"HeatFlu");
447  set_of_solvent_fluxes[it->getMeshsetId()].dAta.data.flag1 = 1;
448  set_of_solvent_fluxes[it->getMeshsetId()].dAta.data.value1 = data[0];
449  //cerr << set_of_solvent_fluxes[it->getMeshsetId()].dAta << endl;
450  rval = m_field.get_moab().get_entities_by_type(
451  it->meshset,MBTRI,set_of_solvent_fluxes[it->getMeshsetId()].tRis,true
452  ); CHKERRQ_MOAB(rval);
454  set_of_solvent_fluxes[it->getMeshsetId()].tRis,MBTRI,"CHEMICAL_LOAD_FLUX_FE"
455  ); CHKERRQ(ierr);
456  }
457  }
458  }
459 
460  //build finite elements
461  ierr = m_field.build_finite_elements(); CHKERRQ(ierr);
462  //build adjacencies
463  ierr = m_field.build_adjacencies(bit_level0); CHKERRQ(ierr);
464  }
465 
466  }
467 
468  // attach tags for each recorder
469  vector<int> tags;
470  tags.push_back(Gel::STRESSTOTAL); // ADOL-C tag used to calculate total stress
471  tags.push_back(Gel::SOLVENTFLUX);
472  tags.push_back(Gel::SOLVENTRATE);
473  tags.push_back(Gel::RESIDUALSTRAINHAT);
474 
475  // Create gel instance and set operators.
476  {
477 
478  gel.constitutiveEquationPtr = boost::shared_ptr<UserGelConstitutiveEquation<adouble> >(
479  new UserGelConstitutiveEquation<adouble>(gel.blockMaterialData)
480  );
481 
482  // Set name of fields which has been choose to approximate spatial
483  // displacements, solvent concentration and internal state variables.
484  Gel::CommonData &common_data = gel.commonData;
485  common_data.spatialPositionName = "SPATIAL_POSITION";
486  common_data.spatialPositionNameDot = "SPATIAL_POSITION_DOT";
487  common_data.muName = "CHEMICAL_LOAD";
488  common_data.muNameDot = "CHEMICAL_LOAD_DOT";
489  common_data.strainHatName = "HAT_EPS";
490  common_data.strainHatNameDot = "HAT_EPS_DOT";
491 
492  // Set operators to calculate field values at integration points, both for
493  // left and right hand side elements.
494  Gel::GelFE *fe_ptr[] = { &gel.feRhs, &gel.feLhs };
495  for(int ss = 0;ss<2;ss++) {
496  fe_ptr[ss]->getOpPtrVector().push_back(
497  new Gel::OpGetDataAtGaussPts("SPATIAL_POSITION",common_data,false,true)
498  );
499  fe_ptr[ss]->getOpPtrVector().push_back(
500  new Gel::OpGetDataAtGaussPts("SPATIAL_POSITION_DOT",common_data,false,true)
501  );
502  fe_ptr[ss]->getOpPtrVector().push_back(
503  new Gel::OpGetDataAtGaussPts("CHEMICAL_LOAD",common_data,true,true)
504  );
505  fe_ptr[ss]->getOpPtrVector().push_back(
506  new Gel::OpGetDataAtGaussPts("HAT_EPS",common_data,true,false,MBTET)
507  );
508  fe_ptr[ss]->getOpPtrVector().push_back(
509  new Gel::OpGetDataAtGaussPts("HAT_EPS_DOT",common_data,true,false,MBTET)
510  );
511  }
512 
513 
514  // Right hand side operators
515  gel.feRhs.getOpPtrVector().push_back(
516  new Gel::OpJacobian(
517  "SPATIAL_POSITION",tags,gel.constitutiveEquationPtr,gel.commonData,true,false
518  )
519  );
520  gel.feRhs.getOpPtrVector().push_back(
521  new Gel::OpRhsStressTotal(gel.commonData)
522  );
523  gel.feRhs.getOpPtrVector().push_back(
524  new Gel::OpRhsSolventFlux(gel.commonData)
525  );
526  gel.feRhs.getOpPtrVector().push_back(
527  new Gel::OpRhsSolventConcetrationDot(gel.commonData)
528  );
529  gel.feRhs.getOpPtrVector().push_back(
530  new Gel::OpRhsStrainHat(gel.commonData)
531  );
532 
533  // Left hand side operators
534  gel.feLhs.getOpPtrVector().push_back(
535  new Gel::OpJacobian(
536  "SPATIAL_POSITION",tags,gel.constitutiveEquationPtr,gel.commonData,false,true
537  )
538  );
539  gel.feLhs.getOpPtrVector().push_back(
540  new Gel::OpLhsdxdx(gel.commonData)
541  );
542  gel.feLhs.getOpPtrVector().push_back(
543  new Gel::OpLhsdxdMu(gel.commonData)
544  );
545  gel.feLhs.getOpPtrVector().push_back(
546  new Gel::OpLhsdxdStrainHat(gel.commonData)
547  );
548  gel.feLhs.getOpPtrVector().push_back(
549  new Gel::OpLhsdStrainHatdStrainHat(gel.commonData)
550  );
551  gel.feLhs.getOpPtrVector().push_back(
552  new Gel::OpLhsdStrainHatdx(gel.commonData)
553  );
554  gel.feLhs.getOpPtrVector().push_back(
555  new Gel::OpLhsdMudMu(gel.commonData)
556  );
557  gel.feLhs.getOpPtrVector().push_back(
558  new Gel::OpLhsdMudx(gel.commonData)
559  );
560 
561  }
562 
563  // Create discrete manager instance
564  DM dm;
565  DMType dm_name = "DMGEL";
566  {
567  ierr = DMRegister_MoFEM(dm_name); CHKERRQ(ierr);
568  ierr = DMCreate(PETSC_COMM_WORLD,&dm);CHKERRQ(ierr);
569  ierr = DMSetType(dm,dm_name);CHKERRQ(ierr);
570  ierr = DMMoFEMCreateMoFEM(dm,&m_field,dm_name,bit_level0); CHKERRQ(ierr);
571  ierr = DMSetFromOptions(dm); CHKERRQ(ierr);
572  //add elements to dm
573  ierr = DMMoFEMAddElement(dm,"GEL_FE"); CHKERRQ(ierr);
574  ierr = DMMoFEMAddElement(dm,"FORCE_FE"); CHKERRQ(ierr);
575  ierr = DMMoFEMAddElement(dm,"PRESSURE_FE"); CHKERRQ(ierr);
576  ierr = DMMoFEMAddElement(dm,"CHEMICAL_LOAD_FLUX_FE"); CHKERRQ(ierr);
577  ierr = DMSetUp(dm); CHKERRQ(ierr);
578  }
579 
580  // Create matrices and vectors used for analysis
581  Vec T,F;
582  Mat A;
583  {
585  ierr = VecDuplicate(T,&F); CHKERRQ(ierr);
587  }
588 
589  // Setting finite element methods for Dirichelt boundary conditions
590  SpatialPositionsBCFEMethodPreAndPostProc spatial_position_bc(
591  m_field,"SPATIAL_POSITION",A,T,F
592  );
593  spatial_position_bc.methodsOp.push_back(new TimeForceScale("-my_displacements_history",false));
595  m_field,"CHEMICAL_LOAD","CHEMICAL_LOAD",A,T,F
596  );
597  concentration_bc.methodsOp.push_back(new TimeForceScale("-my_chemical_load_history",false));
598 
599  // Setting finite element method for applying tractions
600  boost::ptr_map<string,NeummanForcesSurface> neumann_forces;
601  boost::ptr_map<string,NodalForce> nodal_forces;
602  boost::ptr_map<string,EdgeForce> edge_forces;
603  {
604  //forces and pressures on surface
605  ierr = MetaNeummanForces::setMomentumFluxOperators(m_field,neumann_forces,PETSC_NULL,"SPATIAL_POSITION"); CHKERRQ(ierr);
606  //noadl forces
607  ierr = MetaNodalForces::setOperators(m_field,nodal_forces,PETSC_NULL,"SPATIAL_POSITION"); CHKERRQ(ierr);
608  //edge forces
609  ierr = MetaEdgeForces::setOperators(m_field,edge_forces,PETSC_NULL,"SPATIAL_POSITION"); CHKERRQ(ierr);
610  for(
611  boost::ptr_map<string,NeummanForcesSurface>::iterator mit = neumann_forces.begin();
612  mit!=neumann_forces.end();mit++
613  ) {
614  mit->second->methodsOp.push_back(new TimeForceScale("-my_load_history",false));
615  }
616  for(
617  boost::ptr_map<string,NodalForce>::iterator mit = nodal_forces.begin();
618  mit!=nodal_forces.end();mit++
619  ) {
620  mit->second->methodsOp.push_back(new TimeForceScale("-my_load_history",false));
621  }
622  for(
623  boost::ptr_map<string,EdgeForce>::iterator mit = edge_forces.begin();
624  mit!=edge_forces.end();mit++
625  ) {
626  mit->second->methodsOp.push_back(new TimeForceScale("-my_load_history",false));
627  }
628  }
629 
630  // Solvent surface element, flux, convection radiation
631  // TODO: add radiation and convection
632  ThermalElement::MyTriFE solvent_surface_fe(m_field);
633  {
634  map<int,ThermalElement::FluxData>::iterator sit = set_of_solvent_fluxes.begin();
635  for(;sit!=set_of_solvent_fluxes.end();sit++) {
636  // add flux operator
637  solvent_surface_fe.getOpPtrVector().push_back(
638  new ThermalElement::OpHeatFlux("CHEMICAL_LOAD",PETSC_NULL,sit->second,true)
639  );
640  }
641  }
642 
643  // Add finite elements to Time Stepping Solver, using Discrete Manager Interface
644  {
645  //Rhs
646  ierr = DMMoFEMTSSetIFunction(dm,DM_NO_ELEMENT,NULL,&spatial_position_bc,NULL); CHKERRQ(ierr);
647  ierr = DMMoFEMTSSetIFunction(dm,DM_NO_ELEMENT,NULL,&concentration_bc,NULL); CHKERRQ(ierr);
648  ierr = DMMoFEMTSSetIFunction(dm,"GEL_FE",&gel.feRhs,NULL,NULL); CHKERRQ(ierr);
649  ierr = DMMoFEMTSSetIFunction(dm,"CHEMICAL_LOAD_FLUX_FE",&solvent_surface_fe,NULL,NULL); CHKERRQ(ierr);
650  {
651  boost::ptr_map<string,NeummanForcesSurface>::iterator mit = neumann_forces.begin();
652  for(;mit!=neumann_forces.end();mit++) {
653  ierr = DMMoFEMTSSetIFunction(dm,mit->first.c_str(),&mit->second->getLoopFe(),NULL,NULL); CHKERRQ(ierr);
654  }
655  }
656  {
657  boost::ptr_map<string,NodalForce>::iterator fit = nodal_forces.begin();
658  for(;fit!=nodal_forces.end();fit++) {
659  ierr = DMMoFEMTSSetIFunction(dm,fit->first.c_str(),&fit->second->getLoopFe(),NULL,NULL); CHKERRQ(ierr);
660  }
661  }
662  {
663  boost::ptr_map<string,EdgeForce>::iterator fit = edge_forces.begin();
664  for(;fit!=edge_forces.end();fit++) {
665  ierr = DMMoFEMTSSetIFunction(dm,fit->first.c_str(),&fit->second->getLoopFe(),NULL,NULL); CHKERRQ(ierr);
666  }
667  }
668  ierr = DMMoFEMTSSetIFunction(dm,DM_NO_ELEMENT,NULL,NULL,&spatial_position_bc); CHKERRQ(ierr);
669  ierr = DMMoFEMTSSetIFunction(dm,DM_NO_ELEMENT,NULL,NULL,&concentration_bc); CHKERRQ(ierr);
670 
671  //Lhs
672  ierr = DMMoFEMTSSetIJacobian(dm,DM_NO_ELEMENT,NULL,&spatial_position_bc,NULL); CHKERRQ(ierr);
673  ierr = DMMoFEMTSSetIJacobian(dm,DM_NO_ELEMENT,NULL,&concentration_bc,NULL); CHKERRQ(ierr);
674  ierr = DMMoFEMTSSetIJacobian(dm,"GEL_FE",&gel.feLhs,NULL,NULL); CHKERRQ(ierr);
675  ierr = DMMoFEMTSSetIJacobian(dm,DM_NO_ELEMENT,NULL,NULL,&spatial_position_bc); CHKERRQ(ierr);
676  ierr = DMMoFEMTSSetIJacobian(dm,DM_NO_ELEMENT,NULL,NULL,&concentration_bc); CHKERRQ(ierr);
677 
678  }
679 
680 
681  // Create Time Stepping solver
682  TS ts;
683  {
684  ierr = TSCreate(PETSC_COMM_WORLD,&ts); CHKERRQ(ierr);
685  ierr = TSSetType(ts,TSBEULER); CHKERRQ(ierr);
686  }
687 
688  {
689  ierr = DMoFEMMeshToLocalVector(dm,T,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
690  ierr = DMoFEMPreProcessFiniteElements(dm,&spatial_position_bc); CHKERRQ(ierr);
691  ierr = DMoFEMPreProcessFiniteElements(dm,&concentration_bc); CHKERRQ(ierr);
692  ierr = DMoFEMMeshToGlobalVector(dm,T,INSERT_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
693  }
694 
695  // Solve problem
696  {
697 
698  ierr = TSSetIFunction(ts,F,PETSC_NULL,PETSC_NULL); CHKERRQ(ierr);
699  ierr = TSSetIJacobian(ts,A,A,PETSC_NULL,PETSC_NULL); CHKERRQ(ierr);
700 
701  //Monitor
702  Gel::MonitorPostProc post_proc(
703  m_field,"DMGEL","GEL_FE",gel.commonData,gel.constitutiveEquationPtr,tags
704  );
705  TsCtx *ts_ctx;
706  DMMoFEMGetTsCtx(dm,&ts_ctx);
707  {
708  ts_ctx->get_postProcess_to_do_Monitor().push_back(&post_proc);
709  ierr = TSMonitorSet(ts,f_TSMonitorSet,ts_ctx,PETSC_NULL); CHKERRQ(ierr);
710  }
711 
712  double ftime = 1;
713  ierr = TSSetDuration(ts,PETSC_DEFAULT,ftime); CHKERRQ(ierr);
714  ierr = TSSetFromOptions(ts); CHKERRQ(ierr);
715  ierr = TSSetDM(ts,dm); CHKERRQ(ierr);
716  #if PETSC_VERSION_GE(3,7,0)
717  ierr = TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER); CHKERRQ(ierr);
718  #endif
719  ierr = TSSolve(ts,T); CHKERRQ(ierr);
720  ierr = TSGetTime(ts,&ftime); CHKERRQ(ierr);
721  PetscInt steps,snesfails,rejects,nonlinits,linits;
722  ierr = TSGetTimeStepNumber(ts,&steps); CHKERRQ(ierr);
723  ierr = TSGetSNESFailures(ts,&snesfails); CHKERRQ(ierr);
724  ierr = TSGetStepRejections(ts,&rejects); CHKERRQ(ierr);
725  ierr = TSGetSNESIterations(ts,&nonlinits); CHKERRQ(ierr);
726  ierr = TSGetKSPIterations(ts,&linits); CHKERRQ(ierr);
727  PetscPrintf(PETSC_COMM_WORLD,
728  "steps %D (%D rejected, %D SNES fails), ftime %g, nonlinits %D, linits %D\n",
729  steps,rejects,snesfails,ftime,nonlinits,linits
730  );
731  }
732 
733  if(is_atom_test) {
734  double sum = 0;
735  ierr = VecSum(T,&sum); CHKERRQ(ierr);
736  ierr = PetscPrintf(PETSC_COMM_WORLD,"sum = %9.8e\n",sum); CHKERRQ(ierr);
737  double fnorm;
738  ierr = VecNorm(T,NORM_2,&fnorm); CHKERRQ(ierr);
739  ierr = PetscPrintf(PETSC_COMM_WORLD,"fnorm = %9.8e\n",fnorm); CHKERRQ(ierr);
740 
741  if(fabs(sum-7.57417437e+01)>1e-6) {
742  SETERRQ(PETSC_COMM_WORLD,MOFEM_ATOM_TEST_INVALID,"Failed to pass test");
743  }
744  if(fabs(fnorm-3.67577050e+01)>1e-6) {
745  SETERRQ(PETSC_COMM_WORLD,MOFEM_ATOM_TEST_INVALID,"Failed to pass test");
746  }
747  }
748 
749  // Clean and destroy
750  {
751  ierr = DMDestroy(&dm); CHKERRQ(ierr);
752  ierr = VecDestroy(&T); CHKERRQ(ierr);
753  ierr = VecDestroy(&F); CHKERRQ(ierr);
754  ierr = MatDestroy(&A); CHKERRQ(ierr);
755  //ierr = TSDestroy(&ts); CHKERRQ(ierr);
756  }
757 
758  PetscFinalize();
759 
760  return 0;
761 }
#define CHKERRQ_MOAB(a)
check error code of MoAB function
Definition: definitions.h:482
#define DM_NO_ELEMENT
Definition: DMMoFEM.hpp:23
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
Assemble internal force vector.
Definition: Gels.hpp:1154
Force scale operator for reading two columns.
virtual moab::Interface & get_moab()=0
Common data for gel model.
Definition: Gels.hpp:345
User (hackable) Gel model.
BasicMethodsSequence & get_postProcess_to_do_Monitor()
Definition: TsCtx.hpp:87
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMMoFEM.cpp:440
Interface for Time Stepping (TS) solver.
Definition: TsCtx.hpp:27
static char help[]
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
PetscErrorCode DMMoFEMAddElement(DM dm, const char fe_name[])
add element to dm
Definition: DMMMoFEM.cpp:424
define surface element
Implementation of Gel constitutive modelImplementation follows constitutive equation from: ...
Definition: Gels.hpp:74
definition of volume element
Definition: Gels.hpp:379
PetscErrorCode f_TSMonitorSet(TS ts, PetscInt step, PetscReal t, Vec u, void *ctx)
Definition: TsCtx.cpp:187
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:79
Core (interface) class.
Definition: Core.hpp:50
Calculate internal forces for solvent flux.
Definition: Gels.hpp:1208
PetscErrorCode DMCreateMatrix_MoFEM(DM dm, Mat *M)
Definition: DMMMoFEM.cpp:835
Projection of edge entities with one mid-node on hierarchical basis.
virtual MoFEMErrorCode build_finite_elements(int verb=-1)=0
Build finite elementsBuild finite element data structures. Have to be run before problem and adjacenc...
PetscErrorCode DMCreateGlobalVector_MoFEM(DM dm, Vec *g)
DMShellSetCreateGlobalVectorsets the routine to create a global vector associated with the shell DM...
Definition: DMMMoFEM.cpp:815
PetscErrorCode DMMoFEMTSSetIJacobian(DM dm, const std::string &fe_name, boost::shared_ptr< MoFEM::FEMethod > method, boost::shared_ptr< MoFEM::FEMethod > pre_only, boost::shared_ptr< MoFEM::FEMethod > post_only)
set TS Jacobian evaluation function
Definition: DMMMoFEM.cpp:710
Assemble matrix .
Definition: Gels.hpp:1925
virtual MoFEMErrorCode add_field(const std::string &name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_EXCL, int verb=-1)=0
Add field.
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Common.hpp:147
Assemble matrix .
Definition: Gels.hpp:1842
Residual strain hat.
Definition: Gels.hpp:1311
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:136
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Common.hpp:78
CHKERRQ(ierr)
PetscErrorCode DMoFEMMeshToGlobalVector(DM dm, Vec g, InsertMode mode, ScatterMode scatter_mode)
set ghosted vector values on all existing mesh entities
Definition: DMMMoFEM.cpp:451
DEPRECATED typedef DirichletSpatialPositionsBc SpatialPositionsBCFEMethodPreAndPostProc
virtual MoFEMErrorCode add_ents_to_field_by_type(const Range &ents, const EntityType type, const std::string &name, int verb=-1)=0
Add entities to field meshsetThe lower dimension entities are added depending on the space type...
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:95
virtual MoFEMErrorCode set_field_order(const EntityHandle meshset, const EntityType type, const std::string &name, const ApproximationOrder order, int verb=-1)=0
Set order approximation of the entities in the field.
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:140
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:104
Calculating right hand side.
Definition: Gels.hpp:1260
PetscErrorCode DMoFEMPreProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
Definition: DMMMoFEM.cpp:461
Definition: Gels.hpp:2071
DEPRECATED MoFEMErrorCode set_field(const double val, const EntityType type, const std::string &field_name)
use FieldBlas
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: DMMMoFEM.cpp:150
#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...
boost::ptr_vector< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string &name_filed)=0
set finite element field data
PetscErrorCode DMMoFEMTSSetIFunction(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::FEMethod *pre_only, MoFEM::FEMethod *post_only)
set TS implicit function evaluation function
Definition: DMMMoFEM.cpp:671
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:315
virtual MoFEMErrorCode synchronise_entities(Range &ent, int verb=-1)=0
DEPRECATED typedef DirichletSetFieldFromBlock DirichletBCFromBlockSetFEMethodPreAndPostProc
virtual MoFEMErrorCode build_fields(int verb=-1)=0
virtual MoFEMErrorCode loop_dofs(const Problem *problem_ptr, const std::string &field_name, RowColData rc, EntMethod &method, int lower_rank, int upper_rank, int verb=-1)=0
Make a loop over entities.
PetscErrorCode DMMoFEMGetTsCtx(DM dm, MoFEM::TsCtx **ts_ctx)
get MoFEM::TsCtx data structure
Definition: DMMMoFEM.cpp:791
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL)=0
add finite element
static MoFEMErrorCode setOperators(MoFEM::Interface &m_field, boost::ptr_map< std::string, EdgeForce > &edge_forces, Vec F, const std::string field_name)
Set integration point operators.
Definition: EdgeForce.hpp:107
operator for calculate heat flux and assemble to right hand side mofem_thermal_elem ...
virtual bool check_field(const std::string &name) const =0
check if field is in database
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
continuous field
Definition: definitions.h:175
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...
Assemble matrix .
Definition: Gels.hpp:1497
static MoFEMErrorCode addNeumannBCElements(MoFEM::Interface &m_field, const std::string field_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS", Range *intersect_ptr=NULL)
Declare finite element.
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=-1)=0
build adjacencies
static MoFEMErrorCode setMomentumFluxOperators(MoFEM::Interface &m_field, boost::ptr_map< std::string, NeummanForcesSurface > &neumann_forces, Vec F, const std::string field_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS")
Set operators to finite elements calculating right hand side vector.
Assemble matrix .
Definition: Gels.hpp:1402
field with C-1 continuity
Definition: definitions.h:178

Variable Documentation

◆ help

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

Definition at line 72 of file gel_analysis.cpp.