v0.14.0
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>
31 using namespace MoFEM;
32 
33 #include <PostProcOnRefMesh.hpp>
34 #include <Projection10NodeCoordsOnField.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>
43 using namespace GelModule;
44 
45 #include <boost/program_options.hpp>
46 using namespace std;
47 namespace po = boost::program_options;
48 
49 #include <MethodForForceScaling.hpp>
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>
66 namespace bio = boost::iostreams;
67 using bio::tee_device;
68 using bio::stream;
69 
70 
71 
72 static char help[] = "...\n\n";
73 
74 struct BlockData {
75  int oRder;
76  BlockData():
77  oRder(-1) {
78  }
79 };
80 
81 int 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",
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 }
GelModule::Gel::feLhs
GelFE feLhs
Definition: Gels.hpp:436
GelModule::Gel::blockMaterialData
map< int, BlockMaterialData > blockMaterialData
Definition: Gels.hpp:105
MethodForForceScaling.hpp
GelModule::UserGelConstitutiveEquation
User (hackable) Gel model.
Definition: UserGelModel.hpp:39
MYPCOMM_INDEX
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:215
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:454
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
BlockData::oRder
int oRder
Definition: gel_analysis.cpp:75
L2
@ L2
field with C-1 continuity
Definition: definitions.h:88
DirichletBCFromBlockSetFEMethodPreAndPostProc
DEPRECATED typedef DirichletSetFieldFromBlockWithFlags DirichletBCFromBlockSetFEMethodPreAndPostProc
Definition: DirichletBC.hpp:381
NodalForce.hpp
GelModule
Definition: Gels.hpp:27
ThermalElement::OpHeatFlux
operator for calculate heat flux and assemble to right hand side
Definition: ThermalElement.hpp:380
MoFEM.hpp
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:527
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:501
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:1975
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:1171
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:1201
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:857
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
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
GelModule::Gel::feRhs
GelFE feRhs
Definition: Gels.hpp:436
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.
UserGelModel.hpp
Implementation of Gel finite element.
GelModule::Gel::MonitorPostProc
Definition: Gels.hpp:2071
ThermalElement.hpp
Operators and data structures for thermal analysis.
TimeForceScale.hpp
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:539
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:47
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
Projection10NodeCoordsOnField.hpp
MoFEM::DMoFEMPreProcessFiniteElements
PetscErrorCode DMoFEMPreProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
Definition: DMMoFEM.cpp:550
GelModule::Gel::OpGetDataAtGaussPts
Definition: Gels.hpp:444
PostProcOnRefMesh.hpp
Post-process fields on refined mesh.
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:118
GelModule::Gel::OpLhsdMudMu
Assemble matrix .
Definition: Gels.hpp:1842
GelModule::Gel::commonData
CommonData commonData
Definition: Gels.hpp:376
mesh_file_name
char mesh_file_name[255]
Definition: mesh_smoothing.cpp:23
BlockData::BlockData
BlockData()
Definition: gel_analysis.cpp:76
Gels.hpp
Implementation of Gel finite element.
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
main
int main(int argc, char *argv[])
Definition: gel_analysis.cpp:81
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:98
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
_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:1094
std
Definition: enable_if.hpp:5
BLOCKSET
@ BLOCKSET
Definition: definitions.h:148
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:804
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::BlockData
Definition: MeshsetsManager.cpp:752
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:1146
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
EdgeForce.hpp
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::constitutiveEquationPtr
boost::shared_ptr< Gel::ConstitutiveEquation< adouble > > constitutiveEquationPtr
Definition: Gels.hpp:340
DirichletBC.hpp
GelModule::Gel::OpRhsSolventFlux
Calculate internal forces for solvent flux.
Definition: Gels.hpp:1208
SurfacePressure.hpp
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