v0.14.0
MWLS.cpp
Go to the documentation of this file.
1 
2 /** \file MWLS.cpp
3  */
4 
5 /* This file is part of MoFEM.
6  * MoFEM is free software: you can redistribute it and/or modify it under
7  * the terms of the GNU Lesser General Public License as published by the
8  * Free Software Foundation, either version 3 of the License, or (at your
9  * option) any later version.
10  *
11  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
12  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
13  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
14  * License for more details.
15  *
16  * You should have received a copy of the GNU Lesser General Public
17  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
18 
19 #include <MoFEM.hpp>
20 using namespace MoFEM;
21 #include <cholesky.hpp>
22 #include <BasicFiniteElements.hpp>
23 #include <Mortar.hpp>
24 #include <Hooke.hpp>
25 #include <NeoHookean.hpp>
26 #include <CrackFrontElement.hpp>
27 #include <MWLS.hpp>
28 #include <BasicFiniteElements.hpp>
29 #include <GriffithForceElement.hpp>
30 #include <VolumeLengthQuality.hpp>
31 #include <ComplexConstArea.hpp>
32 #include <CrackPropagation.hpp>
33 
34 namespace FractureMechanics {
35 
36 MWLSApprox::MWLSApprox(MoFEM::Interface &m_field, Vec F_lambda,
37  boost::shared_ptr<DofEntity> arc_length_dof)
38  : mField(m_field), mwlsMoab(mwlsCore), F_lambda(PETSC_NULL),
39  arcLengthDof(arc_length_dof), nbBasePolynomials(1), dmFactor(1),
40  maxElemsInDMFactor(2),
41  useGlobalBaseAtMaterialReferenceConfiguration(PETSC_FALSE),
42  useNodalData(PETSC_FALSE), mwlsAnalyticalInternalStressTest(PETSC_FALSE),
43  approxPointCoords(3) {
44 
45  if (!LogManager::checkIfChannelExist("MWLSWorld")) {
46  auto core_log = logging::core::get();
47 
48  core_log->add_sink(
49  LogManager::createSink(LogManager::getStrmWorld(), "MWLSWorld"));
50  core_log->add_sink(
51  LogManager::createSink(LogManager::getStrmSync(), "MWLSSync"));
52  core_log->add_sink(
53  LogManager::createSink(LogManager::getStrmSelf(), "MWLSSelf"));
54 
55  LogManager::setLog("MWLSWorld");
56  LogManager::setLog("MWLSSync");
57  LogManager::setLog("MWLSSelf");
58 
59  MOFEM_LOG_TAG("MWLSWorld", "MWLS");
60  MOFEM_LOG_TAG("MWLSSync", "MWLS");
61  MOFEM_LOG_TAG("MWLSSelf", "MWLS");
62  }
63 
64  MOFEM_LOG("MWLSWorld", Sev::noisy) << "MWLS created";
65 
66  rhoAtGaussPts = boost::make_shared<VectorDouble>();
67  diffRhoAtGaussPts = boost::make_shared<MatrixDouble>();
68  diffDiffRhoAtGaussPts = boost::make_shared<MatrixDouble>();
69 }
70 
71 MoFEMErrorCode MWLSApprox::getMWLSOptions() {
72  MoFEMFunctionBegin;
73  ierr = PetscOptionsBegin(mField.get_comm(), "",
74  "Get options for moving least square approximation",
75  "none");
76  CHKERRQ(ierr);
77  CHKERR PetscOptionsScalar("-mwls_dm",
78  "edge length scaling for influence radius", "",
79  dmFactor, &dmFactor, PETSC_NULL);
80  MOFEM_LOG_C("MWLSWorld", Sev::inform, "### Input parameter: -mwls_dm %6.4e",
81  dmFactor);
82 
83  CHKERR PetscOptionsInt(
84  "-mwls_number_of_base_functions",
85  "1 = constant, 4 = linear, 10 = quadratic approximation", "",
86  nbBasePolynomials, &nbBasePolynomials, PETSC_NULL);
87  MOFEM_LOG_C("MWLSWorld", Sev::inform,
88  "### Input parameter: -mwls_number_of_base_functions %d",
89  nbBasePolynomials);
90 
91  CHKERR PetscOptionsInt("-mwls_max_elems_factor",
92  "Set max elements factor max_nb_elems = "
93  "maxElemsInDMFactor * nbBasePolynomials",
94  "", maxElemsInDMFactor, &maxElemsInDMFactor,
95  PETSC_NULL);
96  MOFEM_LOG_C("MWLSWorld", Sev::inform,
97  "### Input parameter: -mwls_max_elems_factor %d",
98  maxElemsInDMFactor);
99 
100  CHKERR PetscOptionsBool(
101  "-mwls_use_global_base_at_reference_configuration",
102  "true local mwls base functions at reference configuration are used", "",
103  useGlobalBaseAtMaterialReferenceConfiguration,
104  &useGlobalBaseAtMaterialReferenceConfiguration, NULL);
105  CHKERRQ(ierr);
106 
107  CHKERR PetscOptionsBool(
108  "-mwls_use_nodal_data",
109  "true nodal data values are used (without averaging to the midpoint)", "",
110  useNodalData, &useNodalData, NULL);
111  MOFEM_LOG_C("MWLSWorld", Sev::verbose,
112  "### Input parameter: -mwls_use_nodal_data %d", useNodalData);
113 
114  CHKERR PetscOptionsBool(
115  "-mwls_analytical_internal_stress",
116  "use analytical strain field for internal stress mwls test", "",
117  mwlsAnalyticalInternalStressTest, &mwlsAnalyticalInternalStressTest,
118  NULL);
119  MOFEM_LOG_C("MWLSWorld", Sev::verbose,
120  "### Input parameter: -test_mwls_internal_stress_analytical %d",
121  mwlsAnalyticalInternalStressTest);
122 
123  maxThreeDepth = 30;
124  CHKERR PetscOptionsInt("-mwls_tree_depth", "maximal three depths", "",
125  maxThreeDepth, &maxThreeDepth, PETSC_NULL);
126  MOFEM_LOG("MWLSWorld", Sev::verbose)
127  << "Maximal MWLS three depths " << maxThreeDepth;
128  splitsPerDir = 12;
129  CHKERR PetscOptionsInt("-mwls_splits_per_dir", "splits per direction", "",
130  splitsPerDir, &splitsPerDir, PETSC_NULL);
131  MOFEM_LOG("MWLSWorld", Sev::verbose)
132  << "Splits per direction " << splitsPerDir;
133 
134  ierr = PetscOptionsEnd();
135  CHKERRQ(ierr);
136  MoFEMFunctionReturn(0);
137 }
138 
139 MoFEMErrorCode MWLSApprox::loadMWLSMesh(std::string file_name) {
140  MoFEMFunctionBegin;
141 
142  CHKERR getMWLSOptions();
143 
144  const char *option;
145  option = "";
146  CHKERR mwlsMoab.load_file(file_name.c_str(), 0, option);
147  CHKERR mwlsMoab.get_entities_by_dimension(0, 3, levelTets);
148  if (levelTets.empty()) {
149  SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY,
150  "No 3d element in MWLS mesh");
151  }
152 
153  CHKERR mwlsMoab.get_adjacencies(levelTets, 1, true, levelEdges,
154  moab::Interface::UNION);
155 
156  // Create HO node on tetrahedral
157  if (!useNodalData) {
158  EntityHandle meshset;
159  CHKERR mwlsMoab.create_meshset(MESHSET_SET, meshset);
160  CHKERR mwlsMoab.add_entities(meshset, levelTets);
161  CHKERR mwlsMoab.convert_entities(meshset, false, false, true);
162  CHKERR mwlsMoab.delete_entities(&meshset, 1);
163 
164  // Get HO nodes in TET center
165  std::vector<double> edge_length;
166  for (auto tet : levelTets) {
167  int num_nodes;
168  const EntityHandle *conn;
169  CHKERR mwlsMoab.get_connectivity(tet, conn, num_nodes, false);
170  EntityHandle ho_node;
171  EntityType tit_type = mwlsMoab.type_from_handle(tet);
172  CHKERR mwlsMoab.high_order_node(tet, conn, tit_type, ho_node);
173  levelNodes.insert(ho_node);
174  }
175 
176  } else {
177  Range tets;
178  CHKERR mwlsMoab.get_entities_by_dimension(0, 3, tets);
179  CHKERR mwlsMoab.get_connectivity(tets, levelNodes, true);
180  }
181 
182  // Store edge nodes coordinates in FTensor
183  double edge_node_coords[6];
184  FTensor::Tensor1<double *, 3> t_node_edge[2] = {
185  FTensor::Tensor1<double *, 3>(edge_node_coords, &edge_node_coords[1],
186  &edge_node_coords[2]),
187  FTensor::Tensor1<double *, 3>(&edge_node_coords[3], &edge_node_coords[4],
188  &edge_node_coords[5])};
189 
190  FTensor::Index<'i', 3> i;
191 
192  // Get edge lengths
193  maxEdgeL = 0;
194  for (auto edge : levelEdges) {
195  int num_nodes;
196  const EntityHandle *conn;
197  CHKERR mwlsMoab.get_connectivity(edge, conn, num_nodes, true);
198  CHKERR mwlsMoab.get_coords(conn, num_nodes, edge_node_coords);
199  t_node_edge[0](i) -= t_node_edge[1](i);
200  double l = sqrt(t_node_edge[0](i) * t_node_edge[0](i));
201  maxEdgeL = (maxEdgeL < l) ? l : maxEdgeL;
202  }
203 
204  // Create tree and find maximal edge length. Maximal edge length is used
205  // to find neighbours nodes to material point.
206  myTree = boost::make_shared<BVHTree>(&mwlsMoab);
207  {
208  std::ostringstream opts;
209  opts << "MAX_DEPTH=" << maxThreeDepth
210  // << ";MAX_PER_LEAF=" << maxElemsInDMFactor * nbBasePolynomials + 1
211  << ";SPLITS_PER_DIR=" << splitsPerDir;
212  FileOptions tree_opts(opts.str().c_str());
213 
214  MOFEM_LOG_CHANNEL("MWLSWorld");
215  MOFEM_LOG_TAG("MWLSWorld", "MWLS");
216  BOOST_LOG_SCOPED_THREAD_ATTR("Timeline", attrs::timer());
217  MOFEM_LOG("MWLSWorld", Sev::verbose) << "Build tree ... ";
218  CHKERR myTree->build_tree(levelTets, &treeRoot, &tree_opts);
219  MOFEM_LOG("MWLSWorld", Sev::verbose) << "done";
220  }
221 
222  // moab::BoundBox box;
223  // CHKERR myTree->get_bounding_box(box);
224  // cout << box << endl;
225  // myTree->print();
226 
227  auto map_elems_positions = [&]() {
228  std::map<EntityHandle, std::array<double, 3>> elem_coords_map;
229  Range tets;
230  CHKERR mwlsMoab.get_entities_by_dimension(0, 3, tets, false);
231  MatrixDouble elem_coords;
232  for (auto p = tets.pair_begin(); p != tets.pair_end(); ++p) {
233  EntityHandle *conn;
234  int count, verts_per_e;
235  Range slot(p->first, p->second);
236  CHKERR mwlsMoab.connect_iterate(slot.begin(), slot.end(), conn,
237  verts_per_e, count);
238  if (count != (int)slot.size())
239  THROW_MESSAGE("Should be one continuos sequence");
240  elem_coords.resize(verts_per_e, 3, false);
241 
242  for (auto t = 0; t != count; ++t) {
243  CHKERR mwlsMoab.get_coords(conn, verts_per_e,
244  &*elem_coords.data().begin());
245  auto get_center = [&]() {
246  std::array<double, 3> c{0, 0, 0};
247  for (auto d : {0, 1, 2})
248  for (auto n = 0; n < verts_per_e; ++n)
249  c[d] += elem_coords(n, d);
250  for (auto d : {0, 1, 2})
251  c[d] /= verts_per_e;
252  return c;
253  };
254  elem_coords_map.emplace(std::make_pair(p->first + t, get_center()));
255  conn += verts_per_e;
256  }
257  }
258  return elem_coords_map;
259  };
260 
261  elemCentreMap = map_elems_positions();
262 
263  // Set tag for MWLS stress when SSLV116 test is run
264  if (mwlsAnalyticalInternalStressTest) {
265  Tag th;
266  std::array<double, 9> def;
267  def.fill(0);
268  CHKERR mwlsMoab.tag_get_handle("MED_SSLV116", 9, MB_TYPE_DOUBLE, th,
269  MB_TAG_CREAT | MB_TAG_SPARSE, def.data());
270  }
271 
272  MoFEMFunctionReturn(0);
273 }
274 
275 MoFEMErrorCode MWLSApprox::getValuesToNodes(Tag th) {
276  MoFEMFunctionBegin;
277  if (useNodalData)
278  MoFEMFunctionReturnHot(0);
279 
280  CHKERR mwlsMoab.tag_get_handle("MID_NODE", 1, MB_TYPE_HANDLE, thMidNode,
281  MB_TAG_CREAT | MB_TAG_DENSE);
282 
283  int length;
284  CHKERR mwlsMoab.tag_get_length(th, length);
285 
286  VectorDouble vals(length);
287  vals.clear();
288 
289  if (mwlsAnalyticalInternalStressTest == PETSC_FALSE) {
290 
291  string name;
292  CHKERR mwlsMoab.tag_get_name(th, name);
293 
294  std::function<MoFEMErrorCode(EntityHandle)> get_data_to_node;
295  if (name.compare("RHO") == 0) {
296  get_data_to_node = [this, th, &vals](EntityHandle tet) {
297  MoFEMFunctionBeginHot;
298  Range nodes;
299  double avg_val = 0;
300  CHKERR mwlsMoab.get_connectivity(&tet, 1, nodes, true);
301  for (auto &n : nodes) {
302  double tag_data;
303  CHKERR mwlsMoab.tag_get_data(th, &n, 1, &tag_data);
304  avg_val += tag_data;
305  }
306  avg_val /= nodes.size();
307  vals[0] = avg_val;
308  MoFEMFunctionReturnHot(0);
309  };
310  } else {
311  get_data_to_node = [this, th, &vals](EntityHandle tet) {
312  MoFEMFunctionBeginHot;
313  CHKERR mwlsMoab.tag_get_data(th, &tet, 1, &*vals.begin());
314  MoFEMFunctionReturnHot(0);
315  };
316  }
317 
318  // clean tag on nodes
319  CHKERR mwlsMoab.tag_clear_data(th, levelNodes, &*vals.begin());
320  for (auto tet : levelTets) {
321 
322  CHKERR get_data_to_node(tet);
323  int num_nodes;
324  const EntityHandle *conn;
325  CHKERR mwlsMoab.get_connectivity(tet, conn, num_nodes, false);
326  EntityHandle ho_node;
327  EntityType tet_type = mwlsMoab.type_from_handle(tet);
328  CHKERR mwlsMoab.high_order_node(tet, conn, tet_type, ho_node);
329  CHKERR mwlsMoab.tag_set_data(th, &ho_node, 1, &*vals.begin());
330  CHKERR mwlsMoab.tag_set_data(thMidNode, &tet, 1, &ho_node);
331  }
332 
333  } else {
334 
335  PetscBool add_analytical_internal_stress_operators = PETSC_FALSE;
336  CHKERR PetscOptionsGetBool(PETSC_NULL, "",
337  "-add_analytical_internal_stress_operators",
338  &add_analytical_internal_stress_operators, NULL);
339 
340  CHKERR mwlsMoab.tag_clear_data(th, levelNodes, &*vals.begin());
341  for (auto tet : levelTets) {
342  int num_nodes;
343  const EntityHandle *conn;
344  CHKERR mwlsMoab.get_connectivity(tet, conn, num_nodes, false);
345  EntityHandle ho_node;
346  EntityType tet_type = mwlsMoab.type_from_handle(tet);
347  CHKERR mwlsMoab.high_order_node(tet, conn, tet_type, ho_node);
348  CHKERR mwlsMoab.tag_set_data(thMidNode, &tet, 1, &ho_node);
349 
350  if (add_analytical_internal_stress_operators == PETSC_FALSE) {
351  VectorDouble coords(3);
352  CHKERR mwlsMoab.get_coords(&ho_node, 1, &coords[0]);
353  FTensor::Tensor1<FTensor::PackPtr<double *, 1>, 3> t_coords(
354  &coords[0], &coords[1], &coords[2]);
355  auto t_strain = CrackPropagation::analyticalStrainFunction(t_coords);
356 
357  constexpr double K = 166.666666666e9; // K = E / 3 / (1 - 2*nu),
358  // E = 200 GPa, nu = 0.3
359 
360  // Only diagonal, vectorial notation
361  vals[0] = 3 * K * t_strain(0, 0);
362  vals[1] = 3 * K * t_strain(1, 1);
363  vals[2] = 3 * K * t_strain(2, 2);
364  }
365 
366  CHKERR mwlsMoab.tag_set_data(th, &ho_node, 1, &*vals.begin());
367  }
368  }
369 
370  MoFEMFunctionReturn(0);
371 }
372 
373 MoFEMErrorCode MWLSApprox::getInfluenceNodes(double material_coords[3]) {
374  constexpr double eps_overlap = 0.01;
375  const size_t max_nb_elems = maxElemsInDMFactor * nbBasePolynomials + 1;
376  MoFEMFunctionBegin;
377 
378  using DistEntPair = std::pair<double, EntityHandle>;
379  using DistEntMap = multi_index_container<
380  DistEntPair,
381  indexed_by<
382 
383  ordered_non_unique<member<DistEntPair, double, &DistEntPair::first>>
384 
385  >>;
386  DistEntMap dist_map;
387 
388  // FIXME: That is potentisl source of problems, should be sparate
389  // function for setting approximation base point, or simply should be set
390  // outside of that function. Other functions will behave depending how
391  // approxPointCoords is set. And here is somehow hidden from the view,
392  // i.e. easy to overlook this. TODO: create a new function like:
393  // MoFEMErrorCode setApproxBasePt(double material_coords[3]) ?
394 
395  for (int dd : {0, 1, 2})
396  approxPointCoords[dd] = material_coords[dd];
397 
398  if (myTree) {
399  shortenDm = maxEdgeL * dmFactor; // Influence radius
400 
401  auto find_leafs = [&](const auto dm) {
402  leafsVec.clear();
403  leafsDist.clear();
404  CHKERR myTree->distance_search(material_coords, shortenDm, leafsVec,
405  1.0e-10, 1.0e-6, &leafsDist);
406  distLeafs.clear();
407  distLeafs.reserve(leafsDist.size());
408  for (size_t i = 0; i != leafsVec.size(); ++i) {
409  distLeafs.emplace_back(leafsDist[i], leafsVec[i]);
410  }
411  std::sort(
412  distLeafs.begin(), distLeafs.end(),
413  [](const auto &a, const auto &b) { return a.first < b.first; });
414  return distLeafs;
415  };
416 
417  auto get_influence_nodes = [&](const auto &tets) {
418  std::vector<EntityHandle> influence_nodes;
419  if (!useNodalData) {
420  influence_nodes.resize(tets.size());
421  CHKERR mwlsMoab.tag_get_data(thMidNode, &*tets.begin(), tets.size(),
422  &*influence_nodes.begin());
423  } else {
424  Range nodes;
425  CHKERR mwlsMoab.get_connectivity(&*tets.begin(), tets.size(), nodes,
426  true);
427  influence_nodes.clear();
428  influence_nodes.reserve(nodes.size());
429  for (auto n : nodes)
430  influence_nodes.emplace_back(n);
431  }
432  return influence_nodes;
433  };
434 
435  auto to_range = [&](auto ents) {
436  Range r;
437  r.insert_list(ents.begin(), ents.end());
438  return r;
439  };
440 
441  auto save_entities = [&](const std::string name, Range ents) {
442  MoFEMFunctionBeginHot;
443  if (!ents.empty()) {
444  EntityHandle meshset;
445  CHKERR mwlsMoab.create_meshset(MESHSET_SET, meshset);
446  CHKERR mwlsMoab.add_entities(meshset, ents);
447 
448  CHKERR mwlsMoab.write_file(name.c_str(), "VTK", "", &meshset, 1);
449  CHKERR mwlsMoab.delete_entities(&meshset, 1);
450  }
451  MoFEMFunctionReturnHot(0);
452  };
453 
454  auto leafs = find_leafs(shortenDm);
455 
456  treeTets.clear();
457  std::vector<EntityHandle> leafs_tets;
458  moab::BoundBox box;
459 
460  auto get_dm2_from_influence_points = [&]() {
461  FTensor::Index<'i', 3> i;
462  VectorDouble node_coors_vec(3 * influenceNodes.size());
463  CHKERR mwlsMoab.get_coords(&*influenceNodes.begin(),
464  influenceNodes.size(),
465  &*node_coors_vec.begin());
466  FTensor::Tensor1<double, 3> t_approx_point(
467  material_coords[0], material_coords[1], material_coords[2]);
468  FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_node_coords(
469  &node_coors_vec[0], &node_coors_vec[1], &node_coors_vec[2]);
470  double dm2 = 0;
471  for (auto node : influenceNodes) {
472  t_node_coords(i) -= t_approx_point(i);
473  const double r2 = t_node_coords(i) * t_node_coords(i);
474  if (dm2 < r2)
475  dm2 = r2;
476  ++t_node_coords;
477  }
478  return dm2;
479  };
480 
481  const double shorten_dm2 = shortenDm * shortenDm;
482  double set_dm2 = shorten_dm2;
483  auto min_max_dist_it = dist_map.get<0>().begin();
484  FTensor::Index<'i', 3> i;
485 
486  for (auto l : leafs) {
487  if (dist_map.size() < max_nb_elems ||
488  ((l.first * l.first) <
489  (set_dm2 + std::numeric_limits<double>::epsilon()))) {
490 
491  leafsTetsVecLeaf.clear();
492  CHKERR mwlsMoab.get_entities_by_dimension(l.second, 3,
493  leafsTetsVecLeaf, false);
494 
495  leafsTetsCentre.resize(3 * leafsTetsVecLeaf.size());
496  {
497  FTensor::Tensor1<FTensor::PackPtr<double *, 1>, 3> t_c{
498  &leafsTetsCentre[0], &leafsTetsCentre[leafsTetsVecLeaf.size()],
499  &leafsTetsCentre[2 * leafsTetsVecLeaf.size()]};
500 
501  for (auto tet : leafsTetsVecLeaf) {
502  const auto &c = elemCentreMap[tet];
503  for (auto ii : {0, 1, 2}) {
504  t_c(ii) = c[ii] - material_coords[ii];
505  }
506  ++t_c;
507  }
508  }
509 
510  {
511  FTensor::Tensor1<FTensor::PackPtr<double *, 1>, 3> t_c{
512  &leafsTetsCentre[0], &leafsTetsCentre[leafsTetsVecLeaf.size()],
513  &leafsTetsCentre[2 * leafsTetsVecLeaf.size()]};
514 
515  for (auto tet : leafsTetsVecLeaf) {
516 
517  if (dist_map.size() > max_nb_elems)
518  set_dm2 = std::min(min_max_dist_it->first, shorten_dm2);
519 
520  if (t_c(0) < set_dm2 && t_c(1) < set_dm2 && t_c(2) < set_dm2) {
521  double dm2 = t_c(i) * t_c(i);
522  if (dm2 < set_dm2) {
523  min_max_dist_it =
524  dist_map.get<0>().emplace(std::make_pair(dm2, tet)).first;
525  if (dist_map.size() > max_nb_elems) {
526  int nb = std::distance(dist_map.get<0>().begin(),
527  min_max_dist_it);
528  for (; nb < max_nb_elems; ++nb)
529  ++min_max_dist_it;
530  }
531  }
532  }
533 
534  ++t_c;
535  }
536  }
537  } else {
538  break;
539  }
540  }
541 
542  treeTets.clear();
543  treeTets.reserve(dist_map.size());
544  const double dm2 = set_dm2 + std::numeric_limits<double>::epsilon();
545  auto hi_it = dist_map.get<0>().upper_bound(dm2);
546  for (auto it = dist_map.get<0>().begin(); it != hi_it; ++it)
547  treeTets.emplace_back(it->second);
548 
549  if (treeTets.empty()) {
550  for (auto &it : dist_map)
551  MOFEM_LOG("MWLSWorld", Sev::error)
552  << "dm map " << it.first << " " << it.second;
553 
554  MOFEM_LOG("MWLSWorld", Sev::error) << "leafs found " << leafs.size();
555  MOFEM_LOG("MWLSWorld", Sev::error)
556  << "dist map size " << dist_map.size();
557  MOFEM_LOG("MWLSWorld", Sev::error) << "dm2 " << dm2;
558  MOFEM_LOG("MWLSWorld", Sev::error) << "shorten dm is " << shortenDm;
559  MOFEM_LOG_C("MWLSWorld", Sev::error, "point: ( %g %g %g )",
560  approxPointCoords(0), approxPointCoords(1),
561  approxPointCoords(2));
562  SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "No leafs found.");
563  }
564 
565  influenceNodes = get_influence_nodes(treeTets);
566  const auto dm2_from_influence_pts = get_dm2_from_influence_points();
567  shortenDm = (1 + eps_overlap) * sqrt(dm2_from_influence_pts);
568 
569  // CHKERR save_entities("tree_tets.vtk", to_range(influenceNodes));
570  } else {
571  SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "kd-three not build");
572  }
573  MoFEMFunctionReturn(0);
574 }
575 
576 MoFEMErrorCode MWLSApprox::calculateApproxCoeff(bool trim_influence_nodes,
577  bool global_derivatives) {
578  MoFEMFunctionBegin;
579 
580  auto to_range = [&](auto ents) {
581  Range r;
582  r.insert_list(ents.begin(), ents.end());
583  return r;
584  };
585 
586  auto save_entities = [&](const std::string name, Range ents) {
587  MoFEMFunctionBeginHot;
588  if (!ents.empty()) {
589  EntityHandle meshset;
590  CHKERR mwlsMoab.create_meshset(MESHSET_SET, meshset);
591  CHKERR mwlsMoab.add_entities(meshset, ents);
592 
593  CHKERR mwlsMoab.write_file(name.c_str(), "VTK", "", &meshset, 1);
594  CHKERR mwlsMoab.delete_entities(&meshset, 1);
595  }
596  MoFEMFunctionReturnHot(0);
597  };
598 
599  // Determine points in influence volume
600  auto trim_nodes = [&](const double dm) {
601  MoFEMFunctionBeginHot;
602  FTensor::Index<'i', 3> i;
603  FTensor::Tensor1<double, 3> t_approx_point(
604  approxPointCoords[0], approxPointCoords[1], approxPointCoords[2]);
605  const double dm2 = shortenDm * shortenDm;
606  VectorDouble node_coors_vec(3 * influenceNodes.size());
607  CHKERR mwlsMoab.get_coords(&*influenceNodes.begin(), influenceNodes.size(),
608  &*node_coors_vec.begin());
609  FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_node_coords(
610  &node_coors_vec[0], &node_coors_vec[1], &node_coors_vec[2]);
611  decltype(influenceNodes) influence_nodes_tmp;
612  influence_nodes_tmp.reserve(influenceNodes.size());
613  for (auto node : influenceNodes) {
614  t_node_coords(i) -= t_approx_point(i);
615  const double r2 = t_node_coords(i) * t_node_coords(i);
616  if (r2 < dm2)
617  influence_nodes_tmp.emplace_back(node);
618  ++t_node_coords;
619  }
620  influenceNodes.swap(influence_nodes_tmp);
621  MoFEMFunctionReturnHot(0);
622  };
623 
624  auto eval_AB = [this, global_derivatives](const double dm) {
625  MoFEMFunctionBegin;
626  FTensor::Index<'i', 3> i;
627 
628  A.resize(nbBasePolynomials, nbBasePolynomials, false);
629  B.resize(nbBasePolynomials, influenceNodes.size(), false);
630  A.clear();
631  B.clear();
632  for (int d = 0; d != 3; d++) {
633  diffA[d].resize(nbBasePolynomials, nbBasePolynomials, false);
634  diffB[d].resize(nbBasePolynomials, influenceNodes.size(), false);
635  diffA[d].clear();
636  diffB[d].clear();
637  }
638 
639  FTensor::Tensor1<double, 3> t_approx_point(
640  approxPointCoords[0], approxPointCoords[1], approxPointCoords[2]);
641 
642  VectorDouble nodes_coords(3 * influenceNodes.size());
643  CHKERR mwlsMoab.get_coords(&*influenceNodes.begin(), influenceNodes.size(),
644  &*nodes_coords.begin());
645  double *ptr = &*nodes_coords.begin();
646  FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_node_coords(
647  ptr, ptr + 1, ptr + 2);
648 
649  double min_distance = -1;
650  nearestInfluenceNode = 0;
651 
652  for (int ii = 0; ii != influenceNodes.size(); ++ii) {
653  t_node_coords(i) -= t_approx_point(i);
654  t_node_coords(i) /= shortenDm;
655  const double r = sqrt(t_node_coords(i) * t_node_coords(i));
656 
657  // Get the influence node that has minimum distance from the node in the
658  // projected mesh
659  if (!ii || r < min_distance) {
660  min_distance = r;
661  nearestInfluenceNode = influenceNodes[ii];
662  }
663 
664  FTensor::Tensor1<double, 3> t_diff_r;
665  if (r > 0) {
666  t_diff_r(i) = (1. / r) * t_node_coords(i) * (-1 / shortenDm);
667  } else {
668  t_diff_r(i) = 0;
669  }
670  // Weights
671  const double w = evalWeight(r);
672  calculateBase<FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>>(
673  t_node_coords);
674 
675  // iterate base functions
676  // A_kj = sum_ii w * b_k * b_j
677  //
678  // B_k = w * b_k
679  //
680  // diffA_kj = sum_ii diff_w * b_k * b_j + w * diff_b_k * b_j + w * b_k *
681  // diff_b_j
682  //
683  // diffB = diff_w * b_k + w * diff_bk
684  for (int k = 0; k != nbBasePolynomials; k++) {
685  const double wp = w * baseFun[k];
686  for (int j = 0; j <= k; j++) {
687  A(k, j) += wp * baseFun[j];
688  }
689  B(k, ii) = wp;
690  if (global_derivatives) {
691  const double diff_w_r = evalDiffWeight(r);
692  const double diff_wp_r = diff_w_r * baseFun[k];
693  for (int d = 0; d != 3; d++) {
694  const double diff_wp_r_dx = diff_wp_r * t_diff_r(d);
695  for (int j = 0; j <= k; j++) {
696  diffA[d](k, j) += diff_wp_r_dx * baseFun[j];
697  }
698  diffB[d](k, ii) = diff_wp_r_dx;
699  }
700  }
701  }
702  ++t_node_coords;
703  }
704  if (nearestInfluenceNode == 0) {
705  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
706  "Nearest node not found \n");
707  }
708 
709  MoFEMFunctionReturn(0);
710  };
711 
712  // Invert matrix A
713  auto invert_A = [this]() {
714  MoFEMFunctionBegin;
715 
716  auto solve = [&](const size_t nb_base_functions) {
717  A.resize(nb_base_functions, nb_base_functions, true);
718  L.resize(nb_base_functions, nb_base_functions, false);
719  if (cholesky_decompose(A, L)) {
720  return false;
721  } else {
722 
723  MatrixDouble inv_a(nb_base_functions, nb_base_functions, false);
724  for (int k = 0; k != nb_base_functions; k++) {
725  ublas::matrix_column<MatrixDouble> mr(inv_a, k);
726  mr(k) = 1;
727  cholesky_solve(L, mr, ublas::lower());
728  }
729 
730  if (nb_base_functions == nbBasePolynomials) {
731  invA.swap(inv_a);
732  } else {
733  invA.resize(nbBasePolynomials, nbBasePolynomials, false);
734  invA.clear();
735  for (size_t i = 0; i != nb_base_functions; ++i)
736  for (size_t j = 0; j != nb_base_functions; ++j)
737  invA(i, j) = inv_a(i, j);
738  }
739 
740  return true;
741  }
742  };
743 
744  auto throw_error = [&]() {
745  MoFEMFunctionBegin;
746  FTensor::Tensor1<double, 3> t_approx_point(
747  approxPointCoords[0], approxPointCoords[1], approxPointCoords[2]);
748  cerr << "Point: " << t_approx_point << endl;
749  cerr << "Matrix A: " << A << endl;
750  SETERRQ1(mField.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
751  "Failed to invert matrix - solution could be "
752  "to increase dmFactor to bigger value, e.g. -mwls_dm 4, "
753  "or reduce of number of base fuctions "
754  "-mwls_number_of_base_functions 1. The number of nodes "
755  "found in "
756  "radius is %u",
757  influenceNodes.size());
758  MoFEMFunctionReturn(0);
759  };
760 
761  auto solve_one = [&]() {
762  MoFEMFunctionBegin;
763  if (!solve(1))
764  CHKERR throw_error();
765  MoFEMFunctionReturn(0);
766  };
767 
768  testA = A;
769  if (!solve(nbBasePolynomials)) {
770  if (nbBasePolynomials == 10) {
771  if (!solve(4))
772  CHKERR solve_one();
773  } else if (nbBasePolynomials == 4) {
774  CHKERR solve_one();
775  } else {
776  CHKERR throw_error();
777  }
778  }
779 
780  MoFEMFunctionReturn(0);
781  };
782 
783  auto calculate_shape_functions_coefficients = [this]() {
784  MoFEMFunctionBegin;
785  // Calculate base functions (shape functions) coefficients
786  invAB.resize(invA.size1(), B.size2(), false);
787  noalias(invAB) = prod(invA, B);
788  MoFEMFunctionReturn(0);
789  };
790 
791  // CHKERR trim_nodes(shortenDm);
792  // CHKERR save_entities("tree_tets_trim.vtk", to_range(influenceNodes));
793 
794  CHKERR eval_AB(shortenDm);
795  CHKERR invert_A();
796  CHKERR calculate_shape_functions_coefficients();
797 
798  MoFEMFunctionReturn(0);
799 }
800 
801 MoFEMErrorCode MWLSApprox::evaluateApproxFun(double eval_point[3]) {
802  MoFEMFunctionBegin;
803  // Determine points in influence volume
804 
805  auto cal_base_functions_at_eval_point = [this, eval_point](const double dm) {
806  MoFEMFunctionBegin;
807 
808  FTensor::Tensor1<double, 3> t_approx_point(
809  approxPointCoords[0], approxPointCoords[1], approxPointCoords[2]);
810  FTensor::Tensor1<double, 3> t_eval_point(eval_point[0], eval_point[1],
811  eval_point[2]);
812 
813  FTensor::Index<'i', 3> i;
814  FTensor::Tensor1<double, 3> t_delta;
815  t_delta(i) = t_eval_point(i) - t_approx_point(i);
816  t_delta(i) /= dm;
817 
818  CHKERR calculateBase(t_delta);
819  MoFEMFunctionReturn(0);
820  };
821 
822  auto calculate_shape_functions = [this]() {
823  MoFEMFunctionBegin;
824  // Calculate approximation functions at eval point
825  approxFun.resize(influenceNodes.size(), false);
826  noalias(approxFun) = prod(baseFun, invAB);
827  MoFEMFunctionReturn(0);
828  };
829 
830  CHKERR cal_base_functions_at_eval_point(shortenDm);
831  CHKERR calculate_shape_functions();
832 
833  MoFEMFunctionReturn(0);
834 }
835 
836 MoFEMErrorCode MWLSApprox::evaluateFirstDiffApproxFun(double eval_point[3],
837  bool global_derivatives) {
838  MoFEMFunctionBegin;
839 
840  auto cal_diff_base_functions_at_eval_point = [this,
841  eval_point](const double dm) {
842  MoFEMFunctionBegin;
843 
844  FTensor::Tensor1<double, 3> t_approx_point(
845  approxPointCoords[0], approxPointCoords[1], approxPointCoords[2]);
846  FTensor::Tensor1<double, 3> t_eval_point(eval_point[0], eval_point[1],
847  eval_point[2]);
848 
849  FTensor::Index<'i', 3> i;
850  FTensor::Tensor1<double, 3> t_delta;
851  t_delta(i) = t_eval_point(i) - t_approx_point(i);
852  t_delta(i) /= dm;
853 
854  CHKERR calculateDiffBase(t_delta, dm);
855  MoFEMFunctionReturn(0);
856  };
857 
858  auto calculate_global_shape_functions_derivatives = [this]() {
859  MoFEMFunctionBegin;
860 
861  // Calculate derivatives of base functions
862  VectorDouble tmp;
863  tmp.resize(nbBasePolynomials, false);
864  baseFunInvA.resize(invA.size2(), false);
865  noalias(baseFunInvA) = prod(baseFun, invA);
866 
867  for (int d = 0; d != 3; d++) {
868  diffInvA[d].resize(nbBasePolynomials, nbBasePolynomials, false);
869  diffApproxFun.resize(3, influenceNodes.size(), false);
870  }
871 
872  for (int d = 0; d != 3; d++) {
873  MatrixDouble a = prod(diffA[d], invA);
874  noalias(diffInvA[d]) = -prod(invA, a);
875  }
876 
877  // Calculate line derivatives for each coordinate direction
878  for (int d = 0; d != 3; ++d) {
879 
880  ublas::matrix_row<MatrixDouble> mr(diffApproxFun, d);
881 
882  // b,x * A^-1 B
883  noalias(mr) = prod(diffBaseFun[d], invAB);
884 
885  // b * (A^-1 A,x A^-1) * B
886  noalias(tmp) = prod(baseFun, diffInvA[d]);
887  mr += prod(tmp, B);
888 
889  // b * A^-1 * B,x
890  mr += prod(baseFunInvA, diffB[d]);
891  }
892 
893  MoFEMFunctionReturn(0);
894  };
895 
896  auto calculate_local_shape_functions_derivatives = [this]() {
897  MoFEMFunctionBegin;
898 
899  // Calculate derivatives of base functions
900  diffApproxFun.resize(3, influenceNodes.size(), false);
901  // Calculate line derivatives for each coordinate direction
902  for (int d = 0; d != 3; ++d) {
903  ublas::matrix_row<MatrixDouble> mr(diffApproxFun, d);
904  // b,x * A^-1 B
905  noalias(mr) = prod(diffBaseFun[d], invAB);
906  }
907 
908  MoFEMFunctionReturn(0);
909  };
910 
911  CHKERR cal_diff_base_functions_at_eval_point(shortenDm);
912  if (global_derivatives)
913  CHKERR calculate_global_shape_functions_derivatives();
914  else
915  CHKERR calculate_local_shape_functions_derivatives();
916 
917  MoFEMFunctionReturn(0);
918 }
919 
920 MoFEMErrorCode
921 MWLSApprox::evaluateSecondDiffApproxFun(double eval_point[3],
922  bool global_derivatives) {
923  MoFEMFunctionBegin;
924 
925  auto cal_diff_diff_base_functions_at_eval_point =
926  [this, eval_point](const double dm) {
927  MoFEMFunctionBegin;
928 
929  FTensor::Tensor1<double, 3> t_approx_point(
930  approxPointCoords[0], approxPointCoords[1], approxPointCoords[2]);
931  FTensor::Tensor1<double, 3> t_eval_point(eval_point[0], eval_point[1],
932  eval_point[2]);
933 
934  FTensor::Index<'i', 3> i;
935  FTensor::Tensor1<double, 3> t_delta;
936  t_delta(i) = t_eval_point(i) - t_approx_point(i);
937  t_delta(i) /= dm;
938 
939  CHKERR calculateDiffDiffBase(t_delta, dm);
940  MoFEMFunctionReturn(0);
941  };
942 
943  auto calculate_global_shape_functions_second_derivatives = [this]() {
944  MoFEMFunctionBegin;
945  // to be implemented
946  MoFEMFunctionReturn(0);
947  };
948 
949  auto calculate_local_shape_functions_second_derivatives = [this]() {
950  MoFEMFunctionBegin;
951 
952  // Calculate derivatives of base functions
953  // for (int d = 0; d != 3; d++) {
954  diffDiffApproxFun.resize(9, influenceNodes.size(), false);
955  // }
956  // Calculate line derivatives for each coordinate direction
957  for (int n = 0; n != 3; ++n) {
958  for (int d = 0; d != 3; ++d) {
959  const int indx = d + 3 * n;
960  ublas::matrix_row<MatrixDouble> mr(diffDiffApproxFun, indx);
961  // b,xy * A^-1 B
962  noalias(mr) = prod(diffDiffBaseFun[indx], invAB);
963  }
964  }
965  MoFEMFunctionReturn(0);
966  };
967 
968  CHKERR cal_diff_diff_base_functions_at_eval_point(shortenDm);
969  if (global_derivatives) {
970  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented yet");
971  // CHKERR calculate_global_shape_functions_second_derivatives();
972  } else
973  CHKERR calculate_local_shape_functions_second_derivatives();
974 
975  MoFEMFunctionReturn(0);
976 }
977 
978 MoFEMErrorCode MWLSApprox::getTagData(Tag th) {
979  MoFEMFunctionBegin;
980  int length;
981  CHKERR mwlsMoab.tag_get_length(th, length);
982  influenceNodesData.resize(influenceNodes.size(), length, false);
983  CHKERR mwlsMoab.tag_get_data(th, &*influenceNodes.begin(),
984  influenceNodes.size(),
985  &influenceNodesData(0, 0));
986  outData.resize(length, false);
987  outDiffData.resize(3, length, false);
988  outDiffDiffData.resize(9, length, false);
989  MoFEMFunctionReturn(0);
990 }
991 
992 const MWLSApprox::VecVals &MWLSApprox::getDataApprox() {
993  noalias(outData) = prod(approxFun, influenceNodesData);
994  return outData;
995 }
996 
997 const MWLSApprox::MatDiffVals &MWLSApprox::getDiffDataApprox() {
998  noalias(outDiffData) = prod(diffApproxFun, influenceNodesData);
999  return outDiffData;
1000 }
1001 
1002 const MWLSApprox::MatDiffDiffVals &MWLSApprox::getDiffDiffDataApprox() {
1003  noalias(outDiffDiffData) = prod(diffDiffApproxFun, influenceNodesData);
1004  return outDiffDiffData;
1005 }
1006 
1007 MoFEMErrorCode MWLSApprox::getErrorAtNode(
1008  Tag th, double &total_stress_at_node, double &total_stress_error_at_node,
1009  double &hydro_static_error_at_node, double &deviatoric_error_at_node,
1010  double &total_energy, double &total_energy_error,
1011  const double young_modulus, const double poisson_ratio) {
1012  MoFEMFunctionBegin;
1013 
1014  int length;
1015  CHKERR mwlsMoab.tag_get_length(th, length);
1016 
1017  FTensor::Index<'i', 3> i;
1018  FTensor::Index<'j', 3> j;
1019  FTensor::Index<'k', 3> k;
1020  FTensor::Index<'l', 3> l;
1021 
1022  auto get_compliance_matrix = [&]() {
1023  FTensor::Ddg<double, 3, 3> t_C;
1024  t_C(i, j, k, l) = 0;
1025 
1026  const double c = poisson_ratio / young_modulus;
1027  t_C(0, 0, 0, 0) = 1. / young_modulus;
1028  t_C(0, 0, 1, 1) = -c;
1029  t_C(0, 0, 2, 2) = -c;
1030 
1031  t_C(1, 1, 0, 0) = -c;
1032  t_C(1, 1, 1, 1) = 1. / young_modulus;
1033  t_C(1, 1, 2, 2) = -c;
1034 
1035  t_C(2, 2, 0, 0) = -c;
1036  t_C(2, 2, 1, 1) = -c;
1037  t_C(2, 2, 2, 2) = 1. / young_modulus;
1038 
1039  const double d = (1 + poisson_ratio) / young_modulus;
1040  t_C(0, 1, 0, 1) = d;
1041  t_C(0, 2, 0, 2) = d;
1042  t_C(1, 2, 1, 2) = d;
1043 
1044  return t_C;
1045  };
1046 
1047  auto t_C = get_compliance_matrix();
1048 
1049  if (length != outData.size())
1050  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1051  "Tag length is not consistent with outData size. \nMost probably "
1052  "getTagData has not been invoked yet!\n");
1053 
1054  VectorDouble val_at_nearest_influence_node;
1055  val_at_nearest_influence_node.resize(length, false);
1056  val_at_nearest_influence_node.clear();
1057  CHKERR mwlsMoab.tag_get_data(th, &nearestInfluenceNode, 1,
1058  &*val_at_nearest_influence_node.begin());
1059 
1060  // SIXX 0 (00) SIYY 1 (11) SIZZ 2 (22) SIXY 3 (01) SIXZ 4 (02) SIYZ 5 (12)
1061  // T* d00, T* d01, T* d02, T* d11, T* d12, T* d22
1062  auto get_symm_tensor = [](auto &m) {
1063  FTensor::Tensor2_symmetric<double *, 3> t(
1064 
1065  &m(0), &m(3), &m(4),
1066 
1067  &m(1), &m(5),
1068 
1069  &m(2)
1070 
1071  );
1072  return t;
1073  };
1074 
1075  auto t_total_stress_error = FTensor::Tensor2_symmetric<double, 3>();
1076  auto t_stress_at_nearset_node =
1077  get_symm_tensor(val_at_nearest_influence_node);
1078  auto t_mwls_stress = get_symm_tensor(outData);
1079 
1080  auto get_energy_norm = [&](auto &t_s) {
1081  return 0.5 * t_s(i, j) * (t_C(i, j, k, l) * t_s(k, l));
1082  };
1083  auto get_stress_nrm2 = [i, j](auto &t_s) { return t_s(i, j) * t_s(i, j); };
1084 
1085  total_stress_at_node = get_stress_nrm2(t_stress_at_nearset_node);
1086  total_energy = get_energy_norm(t_stress_at_nearset_node);
1087 
1088  t_total_stress_error(i, j) =
1089  t_mwls_stress(i, j) - t_stress_at_nearset_node(i, j);
1090  total_stress_error_at_node = get_stress_nrm2(t_total_stress_error);
1091  total_energy_error = get_energy_norm(t_total_stress_error);
1092 
1093  constexpr auto t_kd = FTensor::Kronecker_Delta_symmetric<int>();
1094  // calculate error of the hydro static part of the stress
1095  double hydro_static_difference =
1096  (t_kd(i, j) * t_total_stress_error(i, j)) / 3.;
1097  hydro_static_error_at_node =
1098  hydro_static_difference * hydro_static_difference;
1099 
1100  // calculate error of the deviator part of the stress
1101  auto t_dev_error = FTensor::Tensor2_symmetric<double, 3>();
1102  t_dev_error(i, j) =
1103  t_total_stress_error(i, j) - hydro_static_difference * t_kd(i, j);
1104  deviatoric_error_at_node = get_stress_nrm2(t_dev_error);
1105 
1106  MoFEMFunctionReturn(0);
1107 }
1108 
1109 template <class T>
1110 MoFEMErrorCode
1111 MWLSApprox::OpMWLSBase<T>::doWork(int side, EntityType type,
1112  DataForcesAndSourcesCore::EntData &data) {
1113  MoFEMFunctionBegin;
1114  if (type != MBVERTEX) {
1115  MoFEMFunctionReturnHot(0);
1116  }
1117 
1118  if (std::is_same<T, VolumeElementForcesAndSourcesCore>::value)
1119  if (mwlsApprox->tetsInBlock.find(this->getFEEntityHandle()) ==
1120  mwlsApprox->tetsInBlock.end()) {
1121  MoFEMFunctionReturnHot(0);
1122  }
1123 
1124  if (std::is_same<T, FaceElementForcesAndSourcesCore>::value)
1125  if (mwlsApprox->trisInBlock.find(this->getFEEntityHandle()) ==
1126  mwlsApprox->trisInBlock.end()) {
1127  MoFEMFunctionReturnHot(0);
1128  }
1129 
1130  const int nb_gauss_pts = data.getN().size1();
1131 
1132  FTensor::Index<'i', 3> i;
1133  FTensor::Index<'j', 3> j;
1134 
1135  mwlsApprox->approxBasePoint.resize(3, nb_gauss_pts, false);
1136  mwlsApprox->approxBasePoint.clear();
1137  mwlsApprox->singularInitialDisplacement.resize(3, nb_gauss_pts, false);
1138  mwlsApprox->singularInitialDisplacement.clear();
1139  mwlsApprox->singularCurrentDisplacement.resize(3, nb_gauss_pts, false);
1140  mwlsApprox->singularCurrentDisplacement.clear();
1141 
1142  if (auto fe_ptr = feSingularPtr.lock()) {
1143 
1144  if (fe_ptr && fe_ptr->singularElement) {
1145 
1146  if (!matPosAtPtsPtr)
1147  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1148  "Matrix for material positions not acclocated");
1149 
1150  auto t_material_positions = getFTensor1FromMat<3>(*matPosAtPtsPtr);
1151  mwlsApprox->mwlsMaterialPositions.resize(matPosAtPtsPtr->size1(),
1152  matPosAtPtsPtr->size2(), false);
1153  auto t_mwls_material_positions =
1154  getFTensor1FromMat<3>(mwlsApprox->mwlsMaterialPositions);
1155 
1156  FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>
1157  t_singular_displacement(&fe_ptr->singularDisp(0, 0),
1158  &fe_ptr->singularDisp(0, 1),
1159  &fe_ptr->singularDisp(0, 2));
1160 
1161  auto t_inital_singular_displacement =
1162  getFTensor1FromMat<3>(mwlsApprox->singularInitialDisplacement);
1163  auto t_current_singular_displacement =
1164  getFTensor1FromMat<3>(mwlsApprox->singularCurrentDisplacement);
1165  if (!matGradPosAtPtsPtr)
1166  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1167  "Matrix for gradient of positions not allocated");
1168  auto t_H = getFTensor2FromMat<3, 3>(*matGradPosAtPtsPtr);
1169  for (int gg = 0; gg != nb_gauss_pts; gg++) {
1170 
1171  t_inital_singular_displacement(i) = t_singular_displacement(i);
1172  t_current_singular_displacement(i) =
1173  t_H(i, j) * t_singular_displacement(j);
1174  t_mwls_material_positions(i) =
1175  t_material_positions(i) + t_current_singular_displacement(i);
1176 
1177  ++t_mwls_material_positions;
1178  ++t_material_positions;
1179  ++t_H;
1180  ++t_singular_displacement;
1181  ++t_inital_singular_displacement;
1182  ++t_current_singular_displacement;
1183  }
1184  } else {
1185  mwlsApprox->mwlsMaterialPositions = *matPosAtPtsPtr;
1186  }
1187  }
1188 
1189  const bool use_global_base =
1190  mwlsApprox->getUseGlobalBaseAtMaterialReferenceConfiguration();
1191 
1192  if (use_global_base) {
1193  auto t_mwls_material_positions =
1194  getFTensor1FromMat<3>(mwlsApprox->mwlsMaterialPositions);
1195  auto t_approx_base_point =
1196  getFTensor1FromMat<3>(mwlsApprox->approxBasePoint);
1197  for (int gg = 0; gg != nb_gauss_pts; gg++) {
1198  t_approx_base_point(i) = t_mwls_material_positions(i);
1199  ++t_approx_base_point;
1200  ++t_mwls_material_positions;
1201  }
1202  } else {
1203 
1204  auto fe_ptr = feSingularPtr.lock();
1205  if (fe_ptr && fe_ptr->singularElement) {
1206  auto t_inital_singular_displacement =
1207  getFTensor1FromMat<3>(mwlsApprox->singularInitialDisplacement);
1208  auto t_approx_base_point =
1209  getFTensor1FromMat<3>(mwlsApprox->approxBasePoint);
1210  auto t_coords_of_gauss_point = this->getFTensor1CoordsAtGaussPts();
1211  for (int gg = 0; gg != nb_gauss_pts; gg++) {
1212  t_approx_base_point(i) =
1213  t_coords_of_gauss_point(i) + t_inital_singular_displacement(i);
1214  ++t_approx_base_point;
1215  ++t_coords_of_gauss_point;
1216  ++t_inital_singular_displacement;
1217  }
1218  } else {
1219  auto t_approx_base_point =
1220  getFTensor1FromMat<3>(mwlsApprox->approxBasePoint);
1221  auto t_coords_of_gauss_point = this->getFTensor1CoordsAtGaussPts();
1222  for (int gg = 0; gg != nb_gauss_pts; gg++) {
1223  t_approx_base_point(i) = t_coords_of_gauss_point(i);
1224  ++t_approx_base_point;
1225  ++t_coords_of_gauss_point;
1226  }
1227  }
1228  }
1229 
1230  CHKERR doMWLSWork(side, type, data);
1231 
1232  MoFEMFunctionReturn(0);
1233 }
1234 
1235 template <class T>
1236 MoFEMErrorCode
1237 MWLSApprox::OpMWLSCalculateBaseCoeffcientsAtGaussPtsTmpl<T>::doMWLSWork(
1238  int side, EntityType type, DataForcesAndSourcesCore::EntData &data) {
1239  MoFEMFunctionBegin;
1240 
1241  FTensor::Index<'i', 3> i;
1242  const int nb_gauss_pts = data.getN().size1();
1243  auto &mwls = this->mwlsApprox;
1244 
1245  if (mwls->invABMap.find(this->getFEEntityHandle()) == mwls->invABMap.end() ||
1246  mwls->influenceNodesMap.find(this->getFEEntityHandle()) ==
1247  mwls->influenceNodesMap.end()) {
1248 
1249  // MOFEM_LOG("MWLSSelf", Sev::noisy) << "Calculate coeffs" << endl;
1250 
1251  const bool use_global_base =
1252  mwls->getUseGlobalBaseAtMaterialReferenceConfiguration();
1253  if (use_global_base)
1254  SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE,
1255  "Coefficients can be precalculated only for local base in "
1256  "reference configuration");
1257 
1258  auto &inv_AB_map = mwls->invABMap[this->getFEEntityHandle()];
1259  inv_AB_map.resize(nb_gauss_pts);
1260  auto &influence_nodes_map =
1261  mwls->influenceNodesMap[this->getFEEntityHandle()];
1262  influence_nodes_map.resize(nb_gauss_pts);
1263  auto &dm_nodes_map = mwls->dmNodesMap[this->getFEEntityHandle()];
1264  dm_nodes_map.resize(nb_gauss_pts);
1265 
1266  auto t_approx_base_point = getFTensor1FromMat<3>(mwls->approxBasePoint);
1267  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
1268 
1269  FTensor::Tensor1<double, 3> t_approx_base_pos;
1270  t_approx_base_pos(i) = t_approx_base_point(i);
1271 
1272  CHKERR mwls->getInfluenceNodes(&t_approx_base_pos(0));
1273  CHKERR mwls->calculateApproxCoeff(true, false);
1274 
1275  influence_nodes_map[gg] = mwls->influenceNodes;
1276  MatrixDouble &invAB = mwls->invAB;
1277  inv_AB_map[gg].resize(invAB.size1(), invAB.size2(), false);
1278  noalias(inv_AB_map[gg]) = invAB;
1279  dm_nodes_map[gg] = mwls->shortenDm;
1280 
1281  ++t_approx_base_point;
1282  }
1283  }
1284 
1285  MoFEMFunctionReturn(0);
1286 }
1287 
1288 template struct MWLSApprox::OpMWLSCalculateBaseCoeffcientsAtGaussPtsTmpl<
1289  VolumeElementForcesAndSourcesCore>;
1290 template struct MWLSApprox::OpMWLSCalculateBaseCoeffcientsAtGaussPtsTmpl<
1291  FaceElementForcesAndSourcesCore>;
1292 
1293 MoFEMErrorCode MWLSApprox::OpMWLSRhoAtGaussUsingPrecalulatedCoeffs::doMWLSWork(
1294  int side, EntityType type, DataForcesAndSourcesCore::EntData &data) {
1295  MoFEMFunctionBegin;
1296 
1297  if (testing)
1298  cerr << "Element " << getFEEntityHandle() << endl;
1299  const bool use_global_base =
1300  mwlsApprox->getUseGlobalBaseAtMaterialReferenceConfiguration();
1301  if (use_global_base)
1302  SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE,
1303  "Coefficients can be precalculated only for local base in "
1304  "reference configuration");
1305 
1306  FTensor::Index<'i', 3> i;
1307 
1308  Tag th;
1309  CHKERR mwlsApprox->mwlsMoab.tag_get_handle(rhoTagName.c_str(), th);
1310 
1311  const int nb_gauss_pts = data.getN().size1();
1312  auto &inv_AB_map = mwlsApprox->invABMap.at(getFEEntityHandle());
1313  auto &influence_nodes_map =
1314  mwlsApprox->influenceNodesMap.at(getFEEntityHandle());
1315  auto &dm_nodes_map = mwlsApprox->dmNodesMap.at(this->getFEEntityHandle());
1316 
1317  VectorDouble &rho = *mwlsApprox->rhoAtGaussPts;
1318  rho.resize(nb_gauss_pts, false);
1319  MatrixDouble &diff_rho = *mwlsApprox->diffRhoAtGaussPts;
1320  diff_rho.resize(3, nb_gauss_pts, false);
1321  MatrixDouble &diff_diff_rho = *mwlsApprox->diffDiffRhoAtGaussPts;
1322  if (calculate2ndDerivative)
1323  diff_diff_rho.resize(9, nb_gauss_pts, false);
1324 
1325  auto t_mwls_material_positions =
1326  getFTensor1FromMat<3>(mwlsApprox->mwlsMaterialPositions);
1327  auto t_approx_base_point = getFTensor1FromMat<3>(mwlsApprox->approxBasePoint);
1328 
1329  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
1330 
1331  FTensor::Tensor1<double, 3> t_mat_pos;
1332  t_mat_pos(i) = t_mwls_material_positions(i);
1333 
1334  if (testing) {
1335  cerr << "material_positions " << t_mwls_material_positions << endl;
1336  }
1337 
1338  for (int dd : {0, 1, 2})
1339  mwlsApprox->approxPointCoords[dd] = t_approx_base_point(dd);
1340 
1341  mwlsApprox->influenceNodes = influence_nodes_map[gg];
1342  mwlsApprox->invAB = inv_AB_map[gg];
1343  mwlsApprox->shortenDm = dm_nodes_map[gg];
1344 
1345  CHKERR mwlsApprox->evaluateApproxFun(&t_mat_pos(0));
1346  CHKERR mwlsApprox->evaluateFirstDiffApproxFun(&t_mat_pos(0),
1347  use_global_base);
1348  if (calculate2ndDerivative)
1349  CHKERR mwlsApprox->evaluateSecondDiffApproxFun(&t_mat_pos(0),
1350  use_global_base);
1351 
1352  CHKERR mwlsApprox->getTagData(th);
1353  const auto &vals = mwlsApprox->getDataApprox();
1354 
1355  rho(gg) = vals[0];
1356  if (calculateDerivative) {
1357  const auto &diff_vals = mwlsApprox->getDiffDataApprox();
1358  for (int ii = 0; ii != 3; ++ii) {
1359  diff_rho(ii, gg) = diff_vals(ii, 0);
1360  }
1361  }
1362 
1363  if (calculate2ndDerivative) {
1364  const auto &diff_diff_vals = mwlsApprox->getDiffDiffDataApprox();
1365  for (int ii = 0; ii != 3; ++ii) {
1366  for (int jj = 0; jj != 3; ++jj)
1367  diff_diff_rho(jj * 3 + ii, gg) = diff_diff_vals(jj * 3 + ii, 0);
1368  }
1369  }
1370 
1371  ++t_mwls_material_positions;
1372  ++t_approx_base_point;
1373  }
1374 
1375  if (testing) {
1376  cerr << "rho " << rho << endl;
1377  cerr << endl;
1378  }
1379 
1380  MoFEMFunctionReturn(0);
1381 }
1382 
1383 MoFEMErrorCode MWLSApprox::OpMWLSRhoAtGaussPts::doMWLSWork(
1384  int side, EntityType type, DataForcesAndSourcesCore::EntData &data) {
1385  MoFEMFunctionBegin;
1386  if (testing)
1387  cerr << "Element " << getFEEntityHandle() << endl;
1388 
1389  FTensor::Index<'i', 3> i;
1390 
1391  const int nb_gauss_pts = data.getN().size1();
1392  const bool use_global_base =
1393  mwlsApprox->getUseGlobalBaseAtMaterialReferenceConfiguration();
1394  Tag th;
1395  CHKERR mwlsApprox->mwlsMoab.tag_get_handle(rhoTagName.c_str(), th);
1396 
1397  VectorDouble &rho = *mwlsApprox->rhoAtGaussPts;
1398  rho.resize(nb_gauss_pts, false);
1399  MatrixDouble &diff_rho = *mwlsApprox->diffRhoAtGaussPts;
1400  diff_rho.resize(3, nb_gauss_pts, false);
1401  MatrixDouble &diff_diff_rho = *mwlsApprox->diffDiffRhoAtGaussPts;
1402  if (calculate2ndDerivative)
1403  diff_diff_rho.resize(9, nb_gauss_pts, false);
1404 
1405  auto t_mwls_material_positions =
1406  getFTensor1FromMat<3>(mwlsApprox->mwlsMaterialPositions);
1407  auto t_approx_base_point = getFTensor1FromMat<3>(mwlsApprox->approxBasePoint);
1408  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
1409 
1410  FTensor::Tensor1<double, 3> t_mat_pos;
1411  t_mat_pos(i) = t_mwls_material_positions(i);
1412  FTensor::Tensor1<double, 3> t_approx_base_pos;
1413  t_approx_base_pos(i) = t_approx_base_point(i);
1414 
1415  if (testing) {
1416  cerr << "material_positions " << t_mwls_material_positions << endl;
1417  cerr << "approx_base_point " << t_approx_base_point << endl;
1418  }
1419  CHKERR mwlsApprox->getInfluenceNodes(&t_approx_base_pos(0));
1420  CHKERR mwlsApprox->calculateApproxCoeff(true, use_global_base);
1421  CHKERR mwlsApprox->evaluateApproxFun(&t_mat_pos(0));
1422  CHKERR mwlsApprox->evaluateFirstDiffApproxFun(&t_mat_pos(0),
1423  use_global_base);
1424  if (calculate2ndDerivative)
1425  CHKERR mwlsApprox->evaluateSecondDiffApproxFun(&t_mat_pos(0),
1426  use_global_base);
1427  CHKERR mwlsApprox->getTagData(th);
1428  const auto &vals = mwlsApprox->getDataApprox();
1429 
1430  rho(gg) = vals[0];
1431  if (calculateDerivative) {
1432  const auto &diff_vals = mwlsApprox->getDiffDataApprox();
1433  for (int ii = 0; ii != 3; ++ii) {
1434  diff_rho(ii, gg) = diff_vals(ii, 0);
1435  }
1436  }
1437 
1438  if (calculate2ndDerivative) {
1439  const auto &diff_diff_vals = mwlsApprox->getDiffDiffDataApprox();
1440  for (int ii = 0; ii != 3; ++ii) {
1441  for (int jj = 0; jj != 3; ++jj)
1442  diff_diff_rho(ii + 3 * jj, gg) = diff_diff_vals(jj, ii);
1443  }
1444  }
1445 
1446  ++t_mwls_material_positions;
1447  ++t_approx_base_point;
1448  }
1449 
1450  if (testing) {
1451  cerr << "rho " << rho << endl;
1452  cerr << endl;
1453  }
1454 
1455  MoFEMFunctionReturn(0);
1456 }
1457 
1458 PetscErrorCode MWLSApprox::OpMWLSRhoPostProcess::doWork(
1459  int side, EntityType type, DataForcesAndSourcesCore::EntData &data) {
1460  MoFEMFunctionBegin;
1461  if (type != MBVERTEX) {
1462  MoFEMFunctionReturnHot(0);
1463  }
1464  auto &mwls_approx = *mwlsApproxPtr;
1465  if (mwls_approx.tetsInBlock.find(getFEEntityHandle()) ==
1466  mwls_approx.tetsInBlock.end()) {
1467  MoFEMFunctionReturnHot(0);
1468  }
1469 
1470  auto &post_proc_mesh = *postProcMeshPtr;
1471  const auto &map_gauss_pts = *mapGaussPtsPtr;
1472  const auto &rho_at_gauss = *mwls_approx.rhoAtGaussPts;
1473  const auto &rho_diff_at_gauss = *mwls_approx.diffRhoAtGaussPts;
1474 
1475  double rho_tag = 0;
1476  VectorDouble diff_rho_tag(3);
1477  diff_rho_tag.clear();
1478  Tag tag_rho;
1479  Tag tag_diff_rho;
1480 
1481  CHKERR post_proc_mesh.tag_get_handle("RHO_APPROX", 1, MB_TYPE_DOUBLE, tag_rho,
1482  MB_TAG_CREAT | MB_TAG_SPARSE, &rho_tag);
1483  CHKERR post_proc_mesh.tag_get_handle(
1484  "RHO_DIFF_APPROX", 3, MB_TYPE_DOUBLE, tag_diff_rho,
1485  MB_TAG_CREAT | MB_TAG_SPARSE, &*diff_rho_tag.begin());
1486 
1487  const int nb_gauss_pts = data.getN().size1();
1488  // for elastic element
1489  // commonData.fieldMap["RHO"].resize(nb_gauss_pts); //delete this
1490 
1491  for (int gg = 0; gg != nb_gauss_pts; gg++) {
1492  const EntityHandle post_proc_node = map_gauss_pts[gg];
1493  rho_tag = rho_at_gauss(gg);
1494  CHKERR post_proc_mesh.tag_set_data(tag_rho, &post_proc_node, 1, &rho_tag);
1495  for (int ii : {0, 1, 2})
1496  diff_rho_tag[ii] = rho_diff_at_gauss(ii, gg);
1497  CHKERR post_proc_mesh.tag_set_data(tag_diff_rho, &post_proc_node, 1,
1498  &*diff_rho_tag.begin());
1499  }
1500  MoFEMFunctionReturn(0);
1501 }
1502 
1503 MoFEMErrorCode
1504 MWLSApprox::OpMWLSStressAtGaussUsingPrecalulatedCoeffs::doMWLSWork(
1505  int side, EntityType type, DataForcesAndSourcesCore::EntData &data) {
1506  MoFEMFunctionBegin;
1507 
1508  if (testing)
1509  cerr << "Element " << getFEEntityHandle() << endl;
1510 
1511  const bool use_global_base =
1512  mwlsApprox->getUseGlobalBaseAtMaterialReferenceConfiguration();
1513  if (use_global_base)
1514  SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE,
1515  "Coefficients can be precalculated only for local base in "
1516  "reference configuration");
1517 
1518  FTensor::Index<'i', 3> i;
1519 
1520  Tag th;
1521  CHKERR mwlsApprox->mwlsMoab.tag_get_handle(stressTagName.c_str(), th);
1522 
1523  const int nb_gauss_pts = data.getN().size1();
1524  auto &inv_AB_map = mwlsApprox->invABMap.at(getFEEntityHandle());
1525  auto &influence_nodes_map =
1526  mwlsApprox->influenceNodesMap.at(getFEEntityHandle());
1527  auto &dm_nodes_map = mwlsApprox->dmNodesMap.at(this->getFEEntityHandle());
1528 
1529  MatrixDouble &stress = mwlsApprox->stressAtGaussPts;
1530  stress.resize(6, nb_gauss_pts, false);
1531  MatrixDouble &diff_stress = mwlsApprox->diffStressAtGaussPts;
1532  diff_stress.resize(18, nb_gauss_pts, false);
1533 
1534  auto t_mwls_material_positions =
1535  getFTensor1FromMat<3>(mwlsApprox->mwlsMaterialPositions);
1536  auto t_approx_base_point = getFTensor1FromMat<3>(mwlsApprox->approxBasePoint);
1537 
1538  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
1539 
1540  FTensor::Tensor1<double, 3> t_mat_pos;
1541  t_mat_pos(i) = t_mwls_material_positions(i);
1542 
1543  if (testing) {
1544  cerr << "material_positions " << t_mwls_material_positions << endl;
1545  }
1546 
1547  for (int dd : {0, 1, 2})
1548  mwlsApprox->approxPointCoords[dd] = t_approx_base_point(dd);
1549 
1550  mwlsApprox->influenceNodes = influence_nodes_map[gg];
1551  mwlsApprox->invAB = inv_AB_map[gg];
1552  mwlsApprox->shortenDm = dm_nodes_map[gg];
1553 
1554  CHKERR mwlsApprox->evaluateApproxFun(&t_mat_pos(0));
1555  CHKERR mwlsApprox->evaluateFirstDiffApproxFun(&t_mat_pos(0),
1556  use_global_base);
1557 
1558  CHKERR mwlsApprox->getTagData(th);
1559  const auto &vals = mwlsApprox->getDataApprox();
1560 
1561  for (int ii = 0; ii != 6; ++ii)
1562  stress(ii, gg) = vals[ii];
1563  if (calculateDerivative) {
1564  const auto &diff_vals = mwlsApprox->getDiffDataApprox();
1565  for (int ii = 0; ii != 6; ++ii) {
1566  for (int jj = 0; jj != 3; ++jj)
1567  diff_stress(jj * 6 + ii, gg) = diff_vals(jj, ii);
1568  }
1569  }
1570  ++t_mwls_material_positions;
1571  ++t_approx_base_point;
1572  }
1573 
1574  if (testing) {
1575  cerr << "stress " << stress << endl;
1576  cerr << endl;
1577  }
1578 
1579  MoFEMFunctionReturn(0);
1580 }
1581 
1582 MoFEMErrorCode MWLSApprox::OpMWLSStressAtGaussPts::doMWLSWork(
1583  int side, EntityType type, DataForcesAndSourcesCore::EntData &data) {
1584  MoFEMFunctionBegin;
1585 
1586  if (testing)
1587  cerr << "Element " << getFEEntityHandle() << endl;
1588 
1589  FTensor::Index<'i', 3> i;
1590 
1591  const int nb_gauss_pts = data.getN().size1();
1592  const bool use_global_base =
1593  mwlsApprox->getUseGlobalBaseAtMaterialReferenceConfiguration();
1594  Tag th;
1595  CHKERR mwlsApprox->mwlsMoab.tag_get_handle(stressTagName.c_str(), th);
1596 
1597  MatrixDouble &stress = mwlsApprox->stressAtGaussPts;
1598  stress.resize(6, nb_gauss_pts, false);
1599  MatrixDouble &diff_stress = mwlsApprox->diffStressAtGaussPts;
1600  diff_stress.resize(18, nb_gauss_pts, false);
1601 
1602  auto t_mwls_material_positions =
1603  getFTensor1FromMat<3>(mwlsApprox->mwlsMaterialPositions);
1604  auto t_approx_base_point = getFTensor1FromMat<3>(mwlsApprox->approxBasePoint);
1605  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
1606 
1607  auto t_material_positions = getFTensor1FromMat<3>(*matPosAtPtsPtr);
1608 
1609  FTensor::Tensor1<double, 3> t_mat_pos;
1610  t_mat_pos(i) = t_mwls_material_positions(i);
1611  FTensor::Tensor1<double, 3> t_approx_base_pos;
1612  t_approx_base_pos(i) = t_approx_base_point(i);
1613 
1614  if (testing) {
1615  cerr << "material_positions " << t_mwls_material_positions << endl;
1616  cerr << "approx_base_point " << t_approx_base_point << endl;
1617  }
1618  CHKERR mwlsApprox->getInfluenceNodes(&t_approx_base_pos(0));
1619  CHKERR mwlsApprox->calculateApproxCoeff(true, use_global_base);
1620  CHKERR mwlsApprox->evaluateApproxFun(&t_mat_pos(0));
1621  CHKERR mwlsApprox->evaluateFirstDiffApproxFun(&t_mat_pos(0),
1622  use_global_base);
1623 
1624  CHKERR mwlsApprox->getTagData(th);
1625  const auto &vals = mwlsApprox->getDataApprox();
1626 
1627  for (int ii = 0; ii != 6; ++ii)
1628  stress(ii, gg) = vals[ii];
1629  if (calculateDerivative) {
1630  const auto &diff_vals = mwlsApprox->getDiffDataApprox();
1631  for (int ii = 0; ii != 6; ++ii) {
1632  for (int jj = 0; jj != 3; ++jj)
1633  diff_stress(jj * 6 + ii, gg) = diff_vals(jj, ii);
1634  }
1635  }
1636  ++t_mwls_material_positions;
1637  ++t_approx_base_point;
1638  }
1639 
1640  if (testing) {
1641  cerr << "stress " << stress << endl;
1642  cerr << endl;
1643  }
1644 
1645  MoFEMFunctionReturn(0);
1646 }
1647 
1648 PetscErrorCode MWLSApprox::OpMWLSStressPostProcess::doWork(
1649  int side, EntityType type, DataForcesAndSourcesCore::EntData &data) {
1650  MoFEMFunctionBegin;
1651  if (type != MBVERTEX) {
1652  MoFEMFunctionReturnHot(0);
1653  }
1654  auto &mwls_approx = *mwlsApproxPtr;
1655  if (mwls_approx.tetsInBlock.find(getFEEntityHandle()) ==
1656  mwls_approx.tetsInBlock.end()) {
1657  MoFEMFunctionReturnHot(0);
1658  }
1659 
1660  auto &post_proc_mesh = *postProcMeshPtr;
1661  const auto &map_gauss_pts = *mapGaussPtsPtr;
1662  const auto &stress = mwls_approx.stressAtGaussPts;
1663 
1664  VectorDouble9 stress_tag(9);
1665  stress_tag.clear();
1666 
1667  Tag tag_stress;
1668  CHKERR post_proc_mesh.tag_get_handle("STRESS_APPROX", 9, MB_TYPE_DOUBLE,
1669  tag_stress, MB_TAG_CREAT | MB_TAG_SPARSE,
1670  &*stress_tag.begin());
1671 
1672  const int nb_gauss_pts = data.getN().size1();
1673  for (int gg = 0; gg != nb_gauss_pts; gg++) {
1674  const EntityHandle post_proc_node = map_gauss_pts[gg];
1675  for (int ii = 0; ii != 6; ++ii)
1676  stress_tag[ii] = stress(ii, gg);
1677  // this can be confusing while postprocessing since paraview by default is
1678  // showing the magnitude of a tensor
1679  CHKERR post_proc_mesh.tag_set_data(tag_stress, &post_proc_node, 1,
1680  &*stress_tag.begin());
1681  }
1682  MoFEMFunctionReturn(0);
1683 }
1684 
1685 MoFEMErrorCode MWLSApprox::OpMWLSStressAndErrorsAtGaussPts::doWork(
1686  int side, EntityType type, DataForcesAndSourcesCore::EntData &data) {
1687  MoFEMFunctionBegin;
1688 
1689  if (type != MBVERTEX)
1690  MoFEMFunctionReturnHot(0);
1691 
1692  FTensor::Index<'i', 3> i;
1693  FTensor::Index<'j', 3> j;
1694 
1695  auto &mwls_approx = *mwlsApproxPtr;
1696  if (mwls_approx.tetsInBlock.find(getFEEntityHandle()) ==
1697  mwls_approx.tetsInBlock.end()) {
1698  MoFEMFunctionReturnHot(0);
1699  }
1700 
1701  const int nb_gauss_pts = data.getN().size1();
1702  // SIXX 0 (00) SIYY 1 (11) SIZZ 2 (22) SIXY 3 (01) SIXZ 4 (02) SIYZ 5 (12)
1703  // T* d00, T* d01, T* d02, T* d11, T* d12, T* d22
1704  MatrixDouble &stress = mwlsApproxPtr->stressAtGaussPts;
1705  FTensor::Tensor2_symmetric<FTensor::PackPtr<double *, 1>, 3> t_stress(
1706  &stress(0, 0), &stress(3, 0), &stress(4, 0), &stress(1, 0), &stress(5, 0),
1707  &stress(2, 0));
1708 
1709  auto get_symm_tensor = [](auto &m) {
1710  FTensor::Tensor2_symmetric<double *, 3> t(
1711 
1712  &m[0], &m[3], &m[4],
1713 
1714  &m[1], &m[5],
1715 
1716  &m[2]
1717 
1718  );
1719  return t;
1720  };
1721 
1722  auto get_stress_nrm2 = [i, j](auto &t_s) { return t_s(i, j) * t_s(i, j); };
1723 
1724  VectorDouble val_at_nearest_influence_node;
1725  val_at_nearest_influence_node.resize(9, false);
1726  val_at_nearest_influence_node.clear();
1727 
1728  Tag th;
1729  CHKERR mwls_approx.mwlsMoab.tag_get_handle(stressTagName.c_str(), th);
1730 
1731  // This is nasty pice of code, but good for tests. It test value in the
1732  // nearest mode. We should use base functions to evaluate stress at
1733  // arbitrary point. However this is good enough when projected mesh and
1734  // approximated mesh is the same. That is used for tests.
1735  CHKERR mwls_approx.mwlsMoab.tag_get_data(
1736  th, &mwls_approx.nearestInfluenceNode, 1,
1737  &*val_at_nearest_influence_node.data().begin());
1738  auto t_stress_at_nearset_node =
1739  get_symm_tensor(val_at_nearest_influence_node);
1740 
1741  auto t_w = getFTensor0IntegrationWeight();
1742  auto vol = getVolume();
1743 
1744  double int_stress = vol * get_stress_nrm2(t_stress_at_nearset_node);
1745  double stress_error = 0;
1746  double hydrostatic_error = 0;
1747  double deviatoric_error = 0;
1748 
1749  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
1750  const double a = t_w * vol;
1751 
1752  FTensor::Tensor2_symmetric<double, 3> t_diff_stress;
1753  t_diff_stress(i, j) = t_stress(i, j) - t_stress_at_nearset_node(i, j);
1754  stress_error += a * get_stress_nrm2(t_diff_stress);
1755 
1756  constexpr auto t_kd = FTensor::Kronecker_Delta_symmetric<int>();
1757  // calculate error of the hydro static part of the stress
1758  const double hydrostatic_error_at_gg =
1759  (t_kd(i, j) * t_diff_stress(i, j)) / 3.;
1760  hydrostatic_error += a * hydrostatic_error_at_gg * hydrostatic_error_at_gg;
1761 
1762  // calculate error of the deviator part of the stress
1763  auto t_dev_error = FTensor::Tensor2_symmetric<double, 3>();
1764  t_dev_error(i, j) = t_diff_stress(i, j) - hydrostatic_error * t_kd(i, j);
1765  deviatoric_error += a * get_stress_nrm2(t_dev_error);
1766 
1767  ++t_w;
1768  ++t_stress;
1769  }
1770 
1771  auto create_tag = [&](auto name) {
1772  constexpr double zero = 0;
1773  Tag tag;
1774  CHKERR mField.get_moab().tag_get_handle(
1775  name, 1, MB_TYPE_DOUBLE, tag, MB_TAG_CREAT | MB_TAG_SPARSE, &zero);
1776  return tag;
1777  };
1778 
1779  auto save_val = [&](auto &&tag, auto &val) {
1780  MoFEMFunctionBegin;
1781  auto const fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
1782  CHKERR mField.get_moab().tag_set_data(tag, &fe_ent, 1, &val);
1783  MoFEMFunctionReturn(0);
1784  };
1785 
1786  stress_error /= int_stress;
1787  hydrostatic_error /= int_stress;
1788  deviatoric_error /= int_stress;
1789  stress_error = sqrt(stress_error);
1790  hydrostatic_error = sqrt(hydrostatic_error);
1791  deviatoric_error = sqrt(deviatoric_error);
1792 
1793  CHKERR save_val(create_tag("INT_STRESS_ERROR"), stress_error);
1794  CHKERR save_val(create_tag("INT_HYDROSTATIC_ERROR"), hydrostatic_error);
1795  CHKERR save_val(create_tag("INT_DEVIATORIC_ERROR"), deviatoric_error);
1796 
1797  MoFEMFunctionReturn(0);
1798 }
1799 
1800 MoFEMErrorCode MWLSApprox::OpMWLSSpatialStressRhs::doWork(
1801  int side, EntityType type, DataForcesAndSourcesCore::EntData &data) {
1802  MoFEMFunctionBegin;
1803  if (mwlsApprox->tetsInBlock.find(getFEEntityHandle()) ==
1804  mwlsApprox->tetsInBlock.end()) {
1805  MoFEMFunctionReturnHot(0);
1806  }
1807  const int nb_dofs = data.getIndices().size();
1808  if (!nb_dofs) {
1809  MoFEMFunctionReturnHot(0);
1810  }
1811  auto t_diff_n = data.getFTensor1DiffN<3>();
1812  // SIXX 0 (00) SIYY 1 (11) SIZZ 2 (22) SIXY 3 (01) SIXZ 4 (02) SIYZ 5 (12)
1813  // T* d00, T* d01, T* d02, T* d11, T* d12, T* d22
1814  MatrixDouble &stress = mwlsApprox->stressAtGaussPts;
1815  FTensor::Tensor2_symmetric<double *, 3> t_stress(
1816  &stress(0, 0), &stress(3, 0), &stress(4, 0), &stress(1, 0), &stress(5, 0),
1817  &stress(2, 0));
1818  FTensor::Index<'i', 3> i;
1819  FTensor::Index<'j', 3> j;
1820  nF.resize(nb_dofs, false);
1821  nF.clear();
1822  const int nb_gauss_pts = data.getN().size1();
1823  auto t_H = getFTensor2FromMat<3, 3>(*matGradPosAtPtsPtr);
1824  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
1825  double det;
1826  CHKERR determinantTensor3by3(t_H, det);
1827  FTensor::Tensor2<double, 3, 3> t_inv_H;
1828  CHKERR invertTensor3by3(t_H, det, t_inv_H);
1829  if (testing) {
1830  cerr << "V*det " << gg << " " << det * getVolume() << endl;
1831  }
1832  double val = getVolume() * getGaussPts()(3, gg) * det;
1833  FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> rhs(&nF[0], &nF[1],
1834  &nF[2]);
1835  int bb = 0;
1836  for (; bb != nb_dofs / 3; ++bb) {
1837  FTensor::Tensor1<double, 3> t_current_diff_n;
1838  t_current_diff_n(j) = t_inv_H(i, j) * t_diff_n(i);
1839  rhs(i) += val * t_stress(i, j) * t_current_diff_n(j);
1840  ++t_diff_n;
1841  ++rhs;
1842  }
1843  for (; bb != data.getN().size2(); ++bb) {
1844  ++t_diff_n;
1845  }
1846  ++t_stress;
1847  ++t_H;
1848  }
1849  Vec f = mwlsApprox->F_lambda != PETSC_NULL ? mwlsApprox->F_lambda
1850  : getFEMethod()->snes_f;
1851  if (mwlsApprox->F_lambda == PETSC_NULL) {
1852  if (auto arc_dof = mwlsApprox->arcLengthDof.lock()) {
1853  nF *= arc_dof->getFieldData();
1854  }
1855  }
1856  if (testing)
1857  cerr << "Element Spatial " << getFEEntityHandle() << " " << nF << endl;
1858  else {
1859  CHKERR VecSetValues(f, nb_dofs, &data.getIndices()[0], &nF[0], ADD_VALUES);
1860  }
1861  MoFEMFunctionReturn(0);
1862 }
1863 
1864 MoFEMErrorCode MWLSApprox::OpMWLSMaterialStressRhs::doWork(
1865  int side, EntityType type, DataForcesAndSourcesCore::EntData &data) {
1866  MoFEMFunctionBegin;
1867  if (mwlsApprox->tetsInBlock.find(getFEEntityHandle()) ==
1868  mwlsApprox->tetsInBlock.end()) {
1869  MoFEMFunctionReturnHot(0);
1870  }
1871  const int nb_dofs = data.getIndices().size();
1872  if (!nb_dofs) {
1873  MoFEMFunctionReturnHot(0);
1874  }
1875  auto t_diff_n = data.getFTensor1DiffN<3>();
1876  // SIXX 0 (00) SIYY 1 (11) SIZZ 2 (22) SIXY 3 (01) SIXZ 4 (02) SIYZ 5 (12)
1877  // T* d00, T* d01, T* d02, T* d11, T* d12, T* d22
1878  MatrixDouble &stress = mwlsApprox->stressAtGaussPts;
1879  FTensor::Tensor2_symmetric<FTensor::PackPtr<double *, 1>, 3> t_stress(
1880  &stress(0, 0), &stress(3, 0), &stress(4, 0), &stress(1, 0), &stress(5, 0),
1881  &stress(2, 0));
1882  FTensor::Index<'i', 3> i;
1883  FTensor::Index<'j', 3> j;
1884  FTensor::Index<'k', 3> k;
1885 
1886  nF.resize(nb_dofs, false);
1887  nF.clear();
1888  const int nb_gauss_pts = data.getN().size1();
1889  auto t_h = getFTensor2FromMat<3, 3>(*spaceGradPosAtPtsPtr);
1890  auto t_H = getFTensor2FromMat<3, 3>(*matGradPosAtPtsPtr);
1891  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
1892  double det;
1893  CHKERR determinantTensor3by3(t_H, det);
1894  FTensor::Tensor2<double, 3, 3> t_inv_H;
1895  CHKERR invertTensor3by3(t_H, det, t_inv_H);
1896  double val = getVolume() * getGaussPts()(3, gg) * det;
1897 
1898  FTensor::Tensor2<double, 3, 3> t_F;
1899  // Calculate gradient of deformation
1900  t_F(i, k) = t_h(i, j) * t_inv_H(j, k);
1901 
1902  constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
1903  FTensor::Tensor2<double, 3, 3> t_eshelby_stress;
1904  t_eshelby_stress(i, k) = -t_F(j, i) * t_stress(j, k);
1905  int bb = 0;
1906  FTensor::Tensor1<double *, 3> rhs(&nF[0], &nF[1], &nF[2], 3);
1907  for (; bb != nb_dofs / 3; ++bb) {
1908  FTensor::Tensor1<double, 3> t_current_diff_n;
1909  t_current_diff_n(j) = t_inv_H(i, j) * t_diff_n(i);
1910  rhs(i) += val * t_eshelby_stress(i, j) * t_current_diff_n(j);
1911  ++t_diff_n;
1912  ++rhs;
1913  }
1914  for (; bb != data.getN().size2(); ++bb) {
1915  ++t_diff_n;
1916  }
1917  ++t_stress;
1918  ++t_H;
1919  ++t_h;
1920  }
1921  int *indices_ptr = &data.getIndices()[0];
1922  if (!forcesOnlyOnEntitiesRow.empty()) {
1923  iNdices.resize(nb_dofs);
1924  noalias(iNdices) = data.getIndices();
1925  indices_ptr = &iNdices[0];
1926  auto &dofs = data.getFieldDofs();
1927  auto dit = dofs.begin();
1928  for (int ii = 0; dit != dofs.end(); dit++, ii++) {
1929  if (auto dof = (*dit)) {
1930  if (forcesOnlyOnEntitiesRow.find(dof->getEnt()) ==
1931  forcesOnlyOnEntitiesRow.end()) {
1932  iNdices[ii] = -1;
1933  }
1934  }
1935  }
1936  }
1937  Vec f = F_lambda != PETSC_NULL ? F_lambda : getFEMethod()->snes_f;
1938  if (mwlsApprox->F_lambda == PETSC_NULL) {
1939  if (auto arc_dof = mwlsApprox->arcLengthDof.lock()) {
1940  nF *= arc_dof->getFieldData();
1941  }
1942  }
1943  CHKERR VecSetOption(f, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
1944  CHKERR VecSetValues(f, nb_dofs, indices_ptr, &nF[0], ADD_VALUES);
1945  MoFEMFunctionReturn(0);
1946 }
1947 
1948 MoFEMErrorCode MWLSApprox::OpMWLSSpatialStressLhs_DX::doWork(
1949  int row_side, int col_side, EntityType row_type, EntityType col_type,
1950  DataForcesAndSourcesCore::EntData &row_data,
1951  DataForcesAndSourcesCore::EntData &col_data) {
1952  MoFEMFunctionBegin;
1953  if (mwlsApprox->tetsInBlock.find(getFEEntityHandle()) ==
1954  mwlsApprox->tetsInBlock.end()) {
1955  MoFEMFunctionReturnHot(0);
1956  }
1957  const int row_nb_dofs = row_data.getIndices().size();
1958  if (!row_nb_dofs) {
1959  MoFEMFunctionReturnHot(0);
1960  }
1961  const int col_nb_dofs = col_data.getIndices().size();
1962  if (!col_nb_dofs) {
1963  MoFEMFunctionReturnHot(0);
1964  }
1965  nA.resize(row_nb_dofs, col_nb_dofs, false);
1966  nA.clear();
1967 
1968  // SIXX 0 (00) SIYY 1 (11) SIZZ 2 (22) SIXY 3 (01) SIXZ 4 (02) SIYZ 5 (12)
1969  MatrixDouble &stress = mwlsApprox->stressAtGaussPts;
1970  FTensor::Tensor2_symmetric<FTensor::PackPtr<double *, 1>, 3> t_stress(
1971  &stress(0, 0), &stress(3, 0), &stress(4, 0), &stress(1, 0), &stress(5, 0),
1972  &stress(2, 0));
1973  MatrixDouble &d_stress = mwlsApprox->diffStressAtGaussPts;
1974  // Note Christod is symmetric on first two indices
1975  FTensor::Christof<FTensor::PackPtr<double *, 1>, 3, 3> t_d_stress(
1976  &d_stress(0 * 6 + 0, 0), &d_stress(0 * 6 + 3, 0), &d_stress(0 * 6 + 4, 0),
1977  &d_stress(0 * 6 + 1, 0), &d_stress(0 * 6 + 5, 0), &d_stress(0 * 6 + 2, 0),
1978 
1979  &d_stress(1 * 6 + 0, 0), &d_stress(1 * 6 + 3, 0), &d_stress(1 * 6 + 4, 0),
1980  &d_stress(1 * 6 + 1, 0), &d_stress(1 * 6 + 5, 0), &d_stress(1 * 6 + 2, 0),
1981 
1982  &d_stress(2 * 6 + 0, 0), &d_stress(2 * 6 + 3, 0), &d_stress(2 * 6 + 4, 0),
1983  &d_stress(2 * 6 + 1, 0), &d_stress(2 * 6 + 5, 0), &d_stress(2 * 6 + 2, 0)
1984 
1985  );
1986 
1987  // scale stress
1988  double lambda = 1;
1989  if (auto arc_dof = mwlsApprox->arcLengthDof.lock()) {
1990  lambda = arc_dof->getFieldData();
1991  }
1992 
1993  auto t_row_diff_n = row_data.getFTensor1DiffN<3>();
1994  const int nb_row_base = row_data.getN().size2();
1995  const int nb_gauss_pts = row_data.getN().size1();
1996  auto t_initial_singular_displacement =
1997  getFTensor1FromMat<3>(mwlsApprox->singularInitialDisplacement);
1998  auto t_H = getFTensor2FromMat<3, 3>(*matGradPosAtPtsPtr);
1999 
2000  FTensor::Index<'i', 3> i;
2001  FTensor::Index<'j', 3> j;
2002  FTensor::Index<'k', 3> k;
2003  FTensor::Index<'m', 3> m;
2004  FTensor::Index<'n', 3> n;
2005 
2006  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
2007 
2008  double det;
2009  CHKERR determinantTensor3by3(t_H, det);
2010  FTensor::Tensor2<double, 3, 3> t_inv_H;
2011  CHKERR invertTensor3by3(t_H, det, t_inv_H);
2012 
2013  double val = getVolume() * getGaussPts()(3, gg) * det;
2014  int rr = 0;
2015  for (; rr != row_nb_dofs / 3; ++rr) {
2016  FTensor::Tensor1<double, 3> t_current_row_diff_n;
2017  t_current_row_diff_n(j) = t_inv_H(i, j) * t_row_diff_n(i);
2018  auto t_col_diff_n = col_data.getFTensor1DiffN<3>(gg, 0);
2019  auto t_col_n = col_data.getFTensor0N(gg, 0);
2020  FTensor::Tensor2<FTensor::PackPtr<double *, 3>, 3, 3> t_mat(
2021  &nA(3 * rr + 0, 0), &nA(3 * rr + 0, 1), &nA(3 * rr + 0, 2),
2022  &nA(3 * rr + 1, 0), &nA(3 * rr + 1, 1), &nA(3 * rr + 1, 2),
2023  &nA(3 * rr + 2, 0), &nA(3 * rr + 2, 1), &nA(3 * rr + 2, 2));
2024  FTensor::Tensor1<double, 3> t_a;
2025  t_a(i) = val * lambda * t_stress(i, j) * t_current_row_diff_n(j);
2026  FTensor::Tensor2<double, 3, 3> t_b;
2027  t_b(i, j) = val * lambda * t_d_stress(i, j, k) * t_current_row_diff_n(k);
2028  FTensor::Tensor3<double, 3, 3, 3> t_c;
2029  t_c(i, m, n) = -val * lambda * (t_stress(i, j) * t_inv_H(n, j)) *
2030  (t_inv_H(k, m) * t_row_diff_n(k));
2031  for (int cc = 0; cc != col_nb_dofs / 3; ++cc) {
2032  t_mat(i, j) += t_a(i) * t_inv_H(k, j) * t_col_diff_n(k);
2033  t_mat(i, j) += t_b(j, i) * t_col_n;
2034  t_mat(i, j) +=
2035  t_b(j, i) * t_col_diff_n(k) * t_initial_singular_displacement(k);
2036  t_mat(i, m) += t_c(i, m, n) * t_col_diff_n(n);
2037  ++t_col_n;
2038  ++t_col_diff_n;
2039  ++t_mat;
2040  }
2041  ++t_row_diff_n;
2042  }
2043  for (; rr != nb_row_base; ++rr) {
2044  ++t_row_diff_n;
2045  }
2046 
2047  // for (int ii = 0; ii != 3; ++ii)
2048  // for (int kk = 0; kk != 3; ++kk)
2049  // for (int ll = kk + 1; ll < 3; ++ll)
2050  // if (t_d_stress(ii, kk, ll) != t_d_stress(ii, ll, kk))
2051  // SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "AAA");
2052 
2053  ++t_H;
2054  ++t_stress;
2055  ++t_d_stress;
2056  ++t_initial_singular_displacement;
2057  }
2058  CHKERR MatSetValues(getFEMethod()->snes_B, row_nb_dofs,
2059  &*row_data.getIndices().begin(), col_nb_dofs,
2060  &*col_data.getIndices().begin(), &*nA.data().begin(),
2061  ADD_VALUES);
2062  MoFEMFunctionReturn(0);
2063 }
2064 
2065 MoFEMErrorCode MWLSApprox::OpMWLSMaterialStressLhs_Dx::doWork(
2066  int row_side, int col_side, EntityType row_type, EntityType col_type,
2067  DataForcesAndSourcesCore::EntData &row_data,
2068  DataForcesAndSourcesCore::EntData &col_data) {
2069  MoFEMFunctionBegin;
2070  if (mwlsApprox->tetsInBlock.find(getFEEntityHandle()) ==
2071  mwlsApprox->tetsInBlock.end()) {
2072  MoFEMFunctionReturnHot(0);
2073  }
2074  const int row_nb_dofs = row_data.getIndices().size();
2075  if (!row_nb_dofs) {
2076  MoFEMFunctionReturnHot(0);
2077  }
2078  const int col_nb_dofs = col_data.getIndices().size();
2079  if (!col_nb_dofs) {
2080  MoFEMFunctionReturnHot(0);
2081  }
2082  nA.resize(row_nb_dofs, col_nb_dofs, false);
2083  nA.clear();
2084  MatrixDouble &stress = mwlsApprox->stressAtGaussPts;
2085  FTensor::Tensor2_symmetric<FTensor::PackPtr<double *, 1>, 3> t_stress(
2086  &stress(0, 0), &stress(3, 0), &stress(4, 0), &stress(1, 0), &stress(5, 0),
2087  &stress(2, 0));
2088 
2089  auto t_row_diff_n = row_data.getFTensor1DiffN<3>();
2090  const int nb_row_base = row_data.getN().size2();
2091  const int nb_gauss_pts = row_data.getN().size1();
2092  auto t_H = getFTensor2FromMat<3, 3>(*matGradPosAtPtsPtr);
2093  auto t_h = getFTensor2FromMat<3, 3>(*spaceGradPosAtPtsPtr);
2094 
2095  FTensor::Index<'i', 3> i;
2096  FTensor::Index<'j', 3> j;
2097  FTensor::Index<'k', 3> k;
2098  FTensor::Index<'l', 3> l;
2099  FTensor::Index<'m', 3> m;
2100  FTensor::Index<'n', 3> n;
2101 
2102  double lambda = 1;
2103  if (auto arc_dof = mwlsApprox->arcLengthDof.lock()) {
2104  lambda = arc_dof->getFieldData();
2105  }
2106 
2107  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
2108 
2109  double det;
2110  CHKERR determinantTensor3by3(t_H, det);
2111  FTensor::Tensor2<double, 3, 3> t_inv_H;
2112  CHKERR invertTensor3by3(t_H, det, t_inv_H);
2113  double val = getVolume() * getGaussPts()(3, gg) * det;
2114 
2115  int rr = 0;
2116  for (; rr != row_nb_dofs / 3; ++rr) {
2117  auto t_col_diff_n = col_data.getFTensor1DiffN<3>(gg, 0);
2118  FTensor::Tensor1<double, 3> t_current_row_diff_n;
2119  t_current_row_diff_n(j) = t_inv_H(i, j) * t_row_diff_n(i);
2120  FTensor::Tensor2<double *, 3, 3> t_mat(
2121  &nA(3 * rr + 0, 0), &nA(3 * rr + 0, 1), &nA(3 * rr + 0, 2),
2122  &nA(3 * rr + 1, 0), &nA(3 * rr + 1, 1), &nA(3 * rr + 1, 2),
2123  &nA(3 * rr + 2, 0), &nA(3 * rr + 2, 1), &nA(3 * rr + 2, 2), 3);
2124  FTensor::Tensor1<double, 3> t_a;
2125  t_a(j) = -val * lambda * t_stress(j, k) * t_current_row_diff_n(k);
2126  int cc = 0;
2127  for (; cc != col_nb_dofs / 3; cc++) {
2128  t_mat(i, j) += t_inv_H(l, i) * t_a(j) * t_col_diff_n(l);
2129  ++t_col_diff_n;
2130  ++t_mat;
2131  }
2132  ++t_row_diff_n;
2133  }
2134  for (; rr != nb_row_base; ++rr) {
2135  ++t_row_diff_n;
2136  }
2137  ++t_stress;
2138  ++t_H;
2139  // ++t_h;
2140  }
2141  int *indices_ptr = &row_data.getIndices()[0];
2142  if (!forcesOnlyOnEntitiesRow.empty()) {
2143  iNdices.resize(row_nb_dofs);
2144  noalias(iNdices) = row_data.getIndices();
2145  indices_ptr = &iNdices[0];
2146  auto &dofs = row_data.getFieldDofs();
2147  auto dit = dofs.begin();
2148  for (int ii = 0; dit != dofs.end(); dit++, ii++) {
2149  if (auto dof = (*dit)) {
2150  if (forcesOnlyOnEntitiesRow.find(dof->getEnt()) ==
2151  forcesOnlyOnEntitiesRow.end()) {
2152  iNdices[ii] = -1;
2153  }
2154  }
2155  }
2156  }
2157  CHKERR MatSetValues(getFEMethod()->snes_B, row_nb_dofs, indices_ptr,
2158  col_nb_dofs, &*col_data.getIndices().begin(),
2159  &*nA.data().begin(), ADD_VALUES);
2160  MoFEMFunctionReturn(0);
2161 }
2162 
2163 MoFEMErrorCode MWLSApprox::OpMWLSMaterialStressLhs_DX::doWork(
2164  int row_side, int col_side, EntityType row_type, EntityType col_type,
2165  DataForcesAndSourcesCore::EntData &row_data,
2166  DataForcesAndSourcesCore::EntData &col_data) {
2167  MoFEMFunctionBegin;
2168  if (mwlsApprox->tetsInBlock.find(getFEEntityHandle()) ==
2169  mwlsApprox->tetsInBlock.end()) {
2170  MoFEMFunctionReturnHot(0);
2171  }
2172  const int row_nb_dofs = row_data.getIndices().size();
2173  if (!row_nb_dofs) {
2174  MoFEMFunctionReturnHot(0);
2175  }
2176  const int col_nb_dofs = col_data.getIndices().size();
2177  if (!col_nb_dofs) {
2178  MoFEMFunctionReturnHot(0);
2179  }
2180  nA.resize(row_nb_dofs, col_nb_dofs, false);
2181  nA.clear();
2182  MatrixDouble &stress = mwlsApprox->stressAtGaussPts;
2183  // SIXX 0 (00) SIYY 1 (11) SIZZ 2 (22) SIXY 3 (01) SIXZ 4 (02) SIYZ 5 (12)
2184  FTensor::Tensor2<FTensor::PackPtr<double *, 1>, 3, 3> t_stress(
2185  &stress(0, 0), &stress(3, 0), &stress(4, 0), &stress(3, 0), &stress(1, 0),
2186  &stress(5, 0), &stress(4, 0), &stress(5, 0), &stress(2, 0));
2187  MatrixDouble &d_stress = mwlsApprox->diffStressAtGaussPts;
2188  FTensor::Tensor3<FTensor::PackPtr<double *, 1>, 3, 3, 3> t_d_stress(
2189  &d_stress(0 * 6 + 0, 0), &d_stress(0 * 6 + 3, 0), &d_stress(0 * 6 + 4, 0),
2190  &d_stress(0 * 6 + 3, 0), &d_stress(0 * 6 + 1, 0), &d_stress(0 * 6 + 5, 0),
2191  &d_stress(0 * 6 + 4, 0), &d_stress(0 * 6 + 5, 0), &d_stress(0 * 6 + 2, 0),
2192 
2193  &d_stress(1 * 6 + 0, 0), &d_stress(1 * 6 + 3, 0), &d_stress(1 * 6 + 4, 0),
2194  &d_stress(1 * 6 + 3, 0), &d_stress(1 * 6 + 1, 0), &d_stress(1 * 6 + 5, 0),
2195  &d_stress(1 * 6 + 4, 0), &d_stress(1 * 6 + 5, 0), &d_stress(1 * 6 + 2, 0),
2196 
2197  &d_stress(2 * 6 + 0, 0), &d_stress(2 * 6 + 3, 0), &d_stress(2 * 6 + 4, 0),
2198  &d_stress(2 * 6 + 3, 0), &d_stress(2 * 6 + 1, 0), &d_stress(2 * 6 + 5, 0),
2199  &d_stress(2 * 6 + 4, 0), &d_stress(2 * 6 + 5, 0), &d_stress(2 * 6 + 2, 0)
2200 
2201  );
2202  auto t_initial_singular_displacement =
2203  getFTensor1FromMat<3>(mwlsApprox->singularInitialDisplacement);
2204 
2205  auto t_row_diff_n = row_data.getFTensor1DiffN<3>();
2206  const int nb_row_base = row_data.getN().size2();
2207  const int nb_gauss_pts = row_data.getN().size1();
2208  auto t_H = getFTensor2FromMat<3, 3>(*matGradPosAtPtsPtr);
2209  auto t_h = getFTensor2FromMat<3, 3>(*spaceGradPosAtPtsPtr);
2210 
2211  FTensor::Index<'i', 3> i;
2212  FTensor::Index<'j', 3> j;
2213  FTensor::Index<'k', 3> k;
2214  FTensor::Index<'l', 3> l;
2215  FTensor::Index<'m', 3> m;
2216  FTensor::Index<'n', 3> n;
2217 
2218  double lambda = 1;
2219  if (auto arc_dof = mwlsApprox->arcLengthDof.lock()) {
2220  lambda = arc_dof->getFieldData();
2221  }
2222 
2223  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
2224 
2225  double det;
2226  CHKERR determinantTensor3by3(t_H, det);
2227  FTensor::Tensor2<double, 3, 3> t_inv_H;
2228  CHKERR invertTensor3by3(t_H, det, t_inv_H);
2229  double val = getVolume() * getGaussPts()(3, gg) * det;
2230 
2231  FTensor::Tensor2<double, 3, 3> t_F;
2232  t_F(i, j) = t_h(i, k) * t_inv_H(k, j);
2233  FTensor::Tensor3<double, 3, 3, 3> t_a;
2234  t_a(i, k, l) = -val * lambda * t_F(j, i) * t_d_stress(l, j, k);
2235  FTensor::Tensor4<double, 3, 3, 3, 3> t_b;
2236  t_b(i, k, m, n) = val * lambda * t_h(j, l) * t_inv_H(l, m) * t_inv_H(n, i) *
2237  t_stress(j, k);
2238  FTensor::Tensor2<double, 3, 3> t_c;
2239  t_c(i, k) = -val * lambda * t_F(j, i) * t_stress(j, k);
2240  int rr = 0;
2241  for (; rr != row_nb_dofs / 3; ++rr) {
2242  FTensor::Tensor1<double, 3> t_current_row_diff_n;
2243  t_current_row_diff_n(j) = t_inv_H(i, j) * t_row_diff_n(i);
2244  auto t_col_diff_n = col_data.getFTensor1DiffN<3>(gg, 0);
2245  auto t_col_n = col_data.getFTensor0N(gg, 0);
2246  FTensor::Tensor2<double *, 3, 3> t_mat(
2247  &nA(3 * rr + 0, 0), &nA(3 * rr + 0, 1), &nA(3 * rr + 0, 2),
2248  &nA(3 * rr + 1, 0), &nA(3 * rr + 1, 1), &nA(3 * rr + 1, 2),
2249  &nA(3 * rr + 2, 0), &nA(3 * rr + 2, 1), &nA(3 * rr + 2, 2), 3);
2250  FTensor::Tensor2<double, 3, 3> t_a_a;
2251  t_a_a(i, l) = t_a(i, k, l) * t_current_row_diff_n(k);
2252  FTensor::Tensor3<double, 3, 3, 3> t_b_b;
2253  t_b_b(i, m, n) = t_b(i, k, m, n) * t_current_row_diff_n(k);
2254  FTensor::Tensor1<double, 3> t_c_c0;
2255  t_c_c0(i) = t_c(i, k) * t_current_row_diff_n(k);
2256  FTensor::Tensor3<double, 3, 3, 3> t_c_c1;
2257  t_c_c1(i, m, n) =
2258  -t_c(i, l) * t_inv_H(k, m) * t_inv_H(n, l) * t_row_diff_n(k);
2259  int cc = 0;
2260  for (; cc != col_nb_dofs / 3; cc++) {
2261  t_mat(i, l) += t_a_a(i, l) * t_col_n;
2262  t_mat(i, l) +=
2263  t_a_a(i, l) * t_col_diff_n(k) * t_initial_singular_displacement(k);
2264  t_mat(i, m) += t_b_b(i, m, n) * t_col_diff_n(n);
2265  t_mat(i, m) += t_c_c0(i) * t_inv_H(n, m) * t_col_diff_n(n);
2266  t_mat(i, m) += t_c_c1(i, m, n) * t_col_diff_n(n);
2267  ++t_col_diff_n;
2268  ++t_col_n;
2269  ++t_mat;
2270  }
2271  ++t_row_diff_n;
2272  }
2273  for (; rr != nb_row_base; ++rr) {
2274  ++t_row_diff_n;
2275  }
2276  ++t_stress;
2277  ++t_d_stress;
2278  ++t_initial_singular_displacement;
2279  ++t_H;
2280  ++t_h;
2281  }
2282  int *indices_ptr = &row_data.getIndices()[0];
2283  if (!forcesOnlyOnEntitiesRow.empty()) {
2284  iNdices.resize(row_nb_dofs);
2285  noalias(iNdices) = row_data.getIndices();
2286  indices_ptr = &iNdices[0];
2287  auto &dofs = row_data.getFieldDofs();
2288  auto dit = dofs.begin();
2289  for (int ii = 0; dit != dofs.end(); dit++, ii++) {
2290  if (auto dof = (*dit)) {
2291  if (forcesOnlyOnEntitiesRow.find(dof->getEnt()) ==
2292  forcesOnlyOnEntitiesRow.end()) {
2293  iNdices[ii] = -1;
2294  }
2295  }
2296  }
2297  }
2298  CHKERR MatSetValues(getFEMethod()->snes_B, row_nb_dofs, indices_ptr,
2299  col_nb_dofs, &*col_data.getIndices().begin(),
2300  &*nA.data().begin(), ADD_VALUES);
2301  MoFEMFunctionReturn(0);
2302 }
2303 } // namespace FractureMechanics
FractureMechanics::MWLSApprox::calculateApproxCoeff
MoFEMErrorCode calculateApproxCoeff(bool trim_influence_nodes, bool global_derivatives)
Calculate approximation function coefficients.
Definition: MWLS.cpp:576
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
MoFEM::K
VectorDouble K
Definition: Projection10NodeCoordsOnField.cpp:125
FractureMechanics::MWLSApprox::evalWeight
double evalWeight(const double r)
Definition: MWLS.hpp:781
FractureMechanics::MWLSApprox::evalDiffWeight
double evalDiffWeight(const double r)
Definition: MWLS.hpp:797
FTensor::Christof
Definition: Christof_value.hpp:9
FractureMechanics::MWLSApprox::OpMWLSRhoAtGaussUsingPrecalulatedCoeffs::testing
const bool testing
Definition: MWLS.hpp:308
FractureMechanics::MWLSApprox::thMidNode
Tag thMidNode
Definition: MWLS.hpp:217
FractureMechanics::MWLSApprox::elemCentreMap
std::map< EntityHandle, std::array< double, 3 > > elemCentreMap
Definition: MWLS.hpp:637
MoFEM::MatSetValues
MoFEMErrorCode MatSetValues(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
Definition: EntitiesFieldData.hpp:1631
FractureMechanics::MWLSApprox::levelNodes
Range levelNodes
Definition: MWLS.hpp:635
FractureMechanics::MWLSApprox::OpMWLSBase< VolumeElementForcesAndSourcesCore >::mwlsApprox
boost::shared_ptr< MWLSApprox > mwlsApprox
Definition: MWLS.hpp:242
PlasticOps::w
double w(double eqiv, double dot_tau, double f, double sigma_y, double sigma_Y)
Definition: PlasticOpsGeneric.hpp:276
FractureMechanics::MWLSApprox::invAB
MatrixDouble invAB
Definition: MWLS.hpp:648
FTensor::Tensor1
Definition: Tensor1_value.hpp:8
MOFEM_LOG_CHANNEL
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
Definition: LogManager.hpp:284
FractureMechanics::MWLSApprox::getDataApprox
const VecVals & getDataApprox()
Definition: MWLS.cpp:992
EntityHandle
FractureMechanics::MWLSApprox::nbBasePolynomials
int nbBasePolynomials
Definition: MWLS.hpp:614
FractureMechanics::MWLSApprox::OpMWLSMaterialStressRhs::doWork
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Definition: MWLS.cpp:1864
young_modulus
double young_modulus
Young modulus.
Definition: plastic.cpp:172
FractureMechanics::MWLSApprox::getDiffDiffDataApprox
const MatDiffDiffVals & getDiffDiffDataApprox()
Definition: MWLS.cpp:1002
FractureMechanics::MWLSApprox::OpMWLSStressAndErrorsAtGaussPts::doWork
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Definition: MWLS.cpp:1685
rho
double rho
Definition: plastic.cpp:191
FractureMechanics::MWLSApprox::OpMWLSSpatialStressLhs_DX::doWork
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
Definition: MWLS.cpp:1948
FractureMechanics::MWLSApprox::diffB
MatrixDouble diffB[3]
Definition: MWLS.hpp:657
FractureMechanics::MWLSApprox::diffDiffApproxFun
MatrixDouble diffDiffApproxFun
Definition: MWLS.hpp:661
FractureMechanics::MWLSApprox::getMWLSOptions
MoFEMErrorCode getMWLSOptions()
get options from command line for MWLS
Definition: MWLS.cpp:71
FractureMechanics::MWLSApprox::F_lambda
Vec F_lambda
Definition: MWLS.hpp:58
MoFEM::CoreInterface::get_comm
virtual MPI_Comm & get_comm() const =0
FractureMechanics::MWLSApprox::getErrorAtNode
MoFEMErrorCode getErrorAtNode(Tag th, double &total_stress_at_node, double &total_stress_error_at_node, double &hydro_static_error_at_node, double &deviatoric_error_at_node, double &total_energy, double &total_energy_error, const double young_modulus, const double poisson_ratio)
Get the norm of the error at nod.
Definition: MWLS.cpp:1007
FractureMechanics::MWLSApprox::B
MatrixDouble B
Definition: MWLS.hpp:645
MWLS.hpp
FractureMechanics::MWLSApprox::myTree
boost::shared_ptr< BVHTree > myTree
Definition: MWLS.hpp:611
Mortar.hpp
FractureMechanics::MWLSApprox::splitsPerDir
int splitsPerDir
Split per direction in the tree.
Definition: MWLS.hpp:626
FractureMechanics::MWLSApprox::loadMWLSMesh
MoFEMErrorCode loadMWLSMesh(std::string file_name)
Definition: MWLS.cpp:139
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
FractureMechanics::MWLSApprox::diffRhoAtGaussPts
boost::shared_ptr< MatrixDouble > diffRhoAtGaussPts
Definition: MWLS.hpp:207
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
FractureMechanics::MWLSApprox::OpMWLSCalculateBaseCoeffcientsAtGaussPtsTmpl::doMWLSWork
MoFEMErrorCode doMWLSWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Do specific work in operator.
Definition: MWLS.cpp:1237
MoFEM::th
Tag th
Definition: Projection10NodeCoordsOnField.cpp:122
FTensor::Kronecker_Delta
Kronecker Delta class.
Definition: Kronecker_Delta.hpp:15
MoFEM.hpp
BasicFiniteElements.hpp
FractureMechanics::MWLSApprox::outData
VecVals outData
Definition: MWLS.hpp:664
FractureMechanics::MWLSApprox::calculateBase
MoFEMErrorCode calculateBase(const T &t_coords)
Definition: MWLS.hpp:675
FractureMechanics::MWLSApprox::evaluateApproxFun
MoFEMErrorCode evaluateApproxFun(double eval_point[3])
evaluate approximation function at material point
Definition: MWLS.cpp:801
FTensor::Tensor2_symmetric
Definition: Tensor2_symmetric_value.hpp:13
FractureMechanics::MWLSApprox::diffApproxFun
MatrixDouble diffApproxFun
Definition: MWLS.hpp:660
MOFEM_IMPOSSIBLE_CASE
@ MOFEM_IMPOSSIBLE_CASE
Definition: definitions.h:35
THROW_MESSAGE
#define THROW_MESSAGE(msg)
Throw MoFEM exception.
Definition: definitions.h:561
FractureMechanics::MWLSApprox::leafsVec
std::vector< EntityHandle > leafsVec
Definition: MWLS.hpp:824
MoFEM::VecSetValues
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
Definition: EntitiesFieldData.hpp:1576
sdf.r
int r
Definition: sdf.py:8
FTensor::Tensor2< double, 3, 3 >
FractureMechanics::MWLSApprox::leafsTetsVecLeaf
std::vector< EntityHandle > leafsTetsVecLeaf
Definition: MWLS.hpp:822
FractureMechanics::MWLSApprox::approxPointCoords
VectorDouble3 approxPointCoords
Definition: MWLS.hpp:214
FractureMechanics::MWLSApprox::outDiffData
MatDiffVals outDiffData
Definition: MWLS.hpp:671
FractureMechanics::MWLSApprox::MatDiffVals
ublas::matrix< double, ublas::row_major, ublas::bounded_array< double, 27 > > MatDiffVals
Definition: MWLS.hpp:178
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
FractureMechanics::MWLSApprox::L
ublas::triangular_matrix< double, ublas::lower > L
Definition: MWLS.hpp:644
FractureMechanics::MWLSApprox::influenceNodesData
MatrixDouble influenceNodesData
Definition: MWLS.hpp:663
FractureMechanics::MWLSApprox::OpMWLSStressAtGaussPts::doMWLSWork
MoFEMErrorCode doMWLSWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Do specific work in operator.
Definition: MWLS.cpp:1582
FractureMechanics::MWLSApprox::influenceNodes
std::vector< EntityHandle > influenceNodes
Definition: MWLS.hpp:636
FractureMechanics::MWLSApprox::approxFun
VectorDouble approxFun
Definition: MWLS.hpp:647
MoFEM::invertTensor3by3
MoFEMErrorCode invertTensor3by3(ublas::matrix< T, L, A > &jac_data, ublas::vector< T, A > &det_data, ublas::matrix< T, L, A > &inv_jac_data)
Calculate inverse of tensor rank 2 at integration points.
Definition: Templates.hpp:1559
FractureMechanics::MWLSApprox::levelEdges
Range levelEdges
Definition: MWLS.hpp:635
c
const double c
speed of light (cm/ns)
Definition: initial_diffusion.cpp:39
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:535
FTensor::Tensor3
Definition: Tensor3_value.hpp:12
FractureMechanics::MWLSApprox::diffInvA
ublas::symmetric_matrix< double, ublas::lower, ublas::row_major, DoubleAllocator > diffInvA[3]
Definition: MWLS.hpp:656
FractureMechanics::MWLSApprox::rhoAtGaussPts
boost::shared_ptr< VectorDouble > rhoAtGaussPts
Definition: MWLS.hpp:206
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
FractureMechanics::MWLSApprox::evaluateSecondDiffApproxFun
MoFEMErrorCode evaluateSecondDiffApproxFun(double eval_point[3], bool global_derivatives)
evaluate 2nd derivative of approximation function at material point
Definition: MWLS.cpp:921
FractureMechanics::MWLSApprox::OpMWLSCalculateBaseCoeffcientsAtGaussPtsTmpl
Definition: MWLS.hpp:277
a
constexpr double a
Definition: approx_sphere.cpp:30
MOFEM_LOG_C
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:311
FractureMechanics::MWLSApprox::A
ublas::symmetric_matrix< double, ublas::lower > A
Definition: MWLS.hpp:641
Hooke.hpp
Implementation of linear elastic material.
convert.type
type
Definition: convert.py:64
FractureMechanics::MWLSApprox::leafsDist
std::vector< double > leafsDist
Definition: MWLS.hpp:825
FractureMechanics::MWLSApprox::OpMWLSStressAtGaussUsingPrecalulatedCoeffs::doMWLSWork
MoFEMErrorCode doMWLSWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Do specific work in operator.
Definition: MWLS.cpp:1504
FractureMechanics::MWLSApprox::dmFactor
double dmFactor
Relative value of weight function radius.
Definition: MWLS.hpp:617
poisson_ratio
double poisson_ratio
Poisson ratio.
Definition: plastic.cpp:173
FractureMechanics::MWLSApprox::maxThreeDepth
int maxThreeDepth
Maximal three depths.
Definition: MWLS.hpp:625
FractureMechanics::MWLSApprox::treeRoot
EntityHandle treeRoot
Definition: MWLS.hpp:632
FractureMechanics::MWLSApprox::OpMWLSRhoPostProcess::doWork
PetscErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Definition: MWLS.cpp:1458
FractureMechanics::MWLSApprox::mwlsMoab
moab::Interface & mwlsMoab
Definition: MWLS.hpp:57
FractureMechanics::MWLSApprox::VecVals
ublas::vector< double, ublas::bounded_array< double, 9 > > VecVals
Definition: MWLS.hpp:175
FractureMechanics::MWLSApprox::testA
ublas::symmetric_matrix< double, ublas::lower > testA
Definition: MWLS.hpp:642
MoFEM::FaceElementForcesAndSourcesCore
Face finite element.
Definition: FaceElementForcesAndSourcesCore.hpp:23
FractureMechanics::MWLSApprox::invA
MatrixDouble invA
Definition: MWLS.hpp:643
FractureMechanics::MWLSApprox::evaluateFirstDiffApproxFun
MoFEMErrorCode evaluateFirstDiffApproxFun(double eval_point[3], bool global_derivatives)
evaluate 1st derivative of approximation function at material point
Definition: MWLS.cpp:836
MOFEM_OPERATION_UNSUCCESSFUL
@ MOFEM_OPERATION_UNSUCCESSFUL
Definition: definitions.h:34
FractureMechanics::MWLSApprox::OpMWLSStressPostProcess::doWork
PetscErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Definition: MWLS.cpp:1648
FractureMechanics::MWLSApprox::OpMWLSMaterialStressLhs_DX::doWork
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
Definition: MWLS.cpp:2163
FractureMechanics::MWLSApprox::calculateDiffDiffBase
MoFEMErrorCode calculateDiffDiffBase(const T &t_coords, const double dm)
Definition: MWLS.hpp:750
FTensor::Tensor4
Definition: Tensor4_value.hpp:18
FractureMechanics::MWLSApprox::outDiffDiffData
MatDiffDiffVals outDiffDiffData
Definition: MWLS.hpp:672
MOFEM_LOG_TAG
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:339
CrackFrontElement.hpp
t
constexpr double t
plate stiffness
Definition: plate.cpp:59
MoFEM::VolumeElementForcesAndSourcesCore
Volume finite element base.
Definition: VolumeElementForcesAndSourcesCore.hpp:26
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
t_kd
constexpr auto t_kd
Definition: free_surface.cpp:137
FractureMechanics::MWLSApprox::OpMWLSBase::doWork
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Definition: MWLS.cpp:1111
FractureMechanics::MWLSApprox::getDiffDataApprox
const MatDiffVals & getDiffDataApprox()
Definition: MWLS.cpp:997
FractureMechanics::MWLSApprox::maxEdgeL
double maxEdgeL
Definition: MWLS.hpp:633
FractureMechanics::CrackPropagation::analyticalStrainFunction
static FTensor::Tensor2_symmetric< double, 3 > analyticalStrainFunction(FTensor::Tensor1< FTensor::PackPtr< double *, 1 >, 3 > &t_coords)
Definition: CrackPropagation.cpp:10030
FractureMechanics::MWLSApprox::baseFunInvA
VectorDouble baseFunInvA
Definition: MWLS.hpp:649
EntData
EntitiesFieldData::EntData EntData
Definition: child_and_parent.cpp:37
FTensor::Index< 'i', 3 >
FractureMechanics::MWLSApprox::baseFun
VectorDouble baseFun
Definition: MWLS.hpp:646
convert.n
n
Definition: convert.py:82
FractureMechanics::MWLSApprox::shortenDm
double shortenDm
Definition: MWLS.hpp:618
FractureMechanics::MWLSApprox::maxElemsInDMFactor
int maxElemsInDMFactor
Definition: MWLS.hpp:620
cholesky.hpp
cholesky decomposition
FractureMechanics::MWLSApprox::mwlsAnalyticalInternalStressTest
PetscBool mwlsAnalyticalInternalStressTest
Definition: MWLS.hpp:630
MoFEM::determinantTensor3by3
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
Definition: Templates.hpp:1511
Range
FTensor::dd
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: ddTensor0.hpp:33
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
FractureMechanics::MWLSApprox::mField
MoFEM::Interface & mField
Definition: MWLS.hpp:55
FractureMechanics::MWLSApprox::OpMWLSSpatialStressRhs::doWork
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Definition: MWLS.cpp:1800
FractureMechanics::MWLSApprox::diffA
ublas::symmetric_matrix< double, ublas::lower, ublas::row_major, DoubleAllocator > diffA[3]
Definition: MWLS.hpp:653
FractureMechanics::MWLSApprox::leafsTetsCentre
std::vector< double > leafsTetsCentre
Definition: MWLS.hpp:823
FractureMechanics::MWLSApprox::diffBaseFun
VectorDouble diffBaseFun[3]
Definition: MWLS.hpp:658
HenckyOps::f
auto f
Definition: HenckyOps.hpp:15
FractureMechanics::MWLSApprox::OpMWLSRhoAtGaussUsingPrecalulatedCoeffs::doMWLSWork
MoFEMErrorCode doMWLSWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Do specific work in operator.
Definition: MWLS.cpp:1293
FractureMechanics::MWLSApprox::useNodalData
PetscBool useNodalData
Definition: MWLS.hpp:629
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
FractureMechanics::MWLSApprox::MatDiffDiffVals
ublas::matrix< double, ublas::row_major, ublas::bounded_array< double, 81 > > MatDiffDiffVals
Definition: MWLS.hpp:181
FTensor::Ddg< double, 3, 3 >
FractureMechanics::MWLSApprox::getInfluenceNodes
std::vector< EntityHandle > & getInfluenceNodes()
Definition: MWLS.hpp:605
MoFEM::Exceptions::ierr
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
FractureMechanics::MWLSApprox::calculateDiffBase
MoFEMErrorCode calculateDiffBase(const T &t_coords, const double dm)
Definition: MWLS.hpp:701
lambda
static double lambda
Definition: incompressible_elasticity.cpp:199
EigenMatrix::Vec
const FTensor::Tensor2< T, Dim, Dim > Vec
Definition: MatrixFunction.hpp:66
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
cholesky_decompose
size_t cholesky_decompose(const MATRIX &A, TRIA &L)
decompose the symmetric positive definit matrix A into product L L^T.
Definition: cholesky.hpp:52
ComplexConstArea.hpp
cholesky_solve
void cholesky_solve(const TRIA &L, VEC &x, ublas::lower)
solve system L L^T x = b inplace
Definition: cholesky.hpp:221
FractureMechanics::MWLSApprox::getTagData
MoFEMErrorCode getTagData(Tag th)
Definition: MWLS.cpp:978
sdf_hertz_2d_axisymm_plane.d
float d
Definition: sdf_hertz_2d_axisymm_plane.py:4
MoFEM::Types::VectorDouble
UBlasVector< double > VectorDouble
Definition: Types.hpp:68
VolumeLengthQuality.hpp
Implementation of Volume-Length-Quality measure with barrier.
m
FTensor::Index< 'm', 3 > m
Definition: shallow_wave.cpp:80
CrackPropagation.hpp
Main class for crack propagation.
FTensor::Kronecker_Delta_symmetric
Kronecker Delta class symmetric.
Definition: Kronecker_Delta.hpp:49
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
FractureMechanics::MWLSApprox::useGlobalBaseAtMaterialReferenceConfiguration
PetscBool useGlobalBaseAtMaterialReferenceConfiguration
Definition: MWLS.hpp:628
FractureMechanics::MWLSApprox::OpMWLSRhoAtGaussUsingPrecalulatedCoeffs::calculateDerivative
const bool calculateDerivative
Definition: MWLS.hpp:306
NeoHookean.hpp
Implementation of Neo-Hookean elastic material.
FractureMechanics::MWLSApprox::levelTets
Range levelTets
Definition: MWLS.hpp:635
FractureMechanics::MWLSApprox::treeTets
std::vector< EntityHandle > treeTets
Definition: MWLS.hpp:821
MoFEM::Types::VectorDouble9
VectorBoundedArray< double, 9 > VectorDouble9
Definition: Types.hpp:96
FractureMechanics::MWLSApprox::OpMWLSRhoAtGaussUsingPrecalulatedCoeffs::rhoTagName
std::string rhoTagName
Definition: MWLS.hpp:305
k
FTensor::Index< 'k', 3 > k
Definition: matrix_function.cpp:20
GriffithForceElement.hpp
FractureMechanics::MWLSApprox::nearestInfluenceNode
EntityHandle nearestInfluenceNode
Definition: MWLS.hpp:639
FractureMechanics::MWLSApprox::getValuesToNodes
MoFEMErrorCode getValuesToNodes(Tag th)
get values form elements to nodes
Definition: MWLS.cpp:275
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
MOFEM_NOT_IMPLEMENTED
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:32
FractureMechanics::MWLSApprox::diffDiffRhoAtGaussPts
boost::shared_ptr< MatrixDouble > diffDiffRhoAtGaussPts
Definition: MWLS.hpp:208
FractureMechanics::MWLSApprox::OpMWLSRhoAtGaussPts::doMWLSWork
MoFEMErrorCode doMWLSWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Do specific work in operator.
Definition: MWLS.cpp:1383
FractureMechanics::MWLSApprox::diffDiffBaseFun
VectorDouble diffDiffBaseFun[9]
Definition: MWLS.hpp:659
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
l
FTensor::Index< 'l', 3 > l
Definition: matrix_function.cpp:21
FractureMechanics::MWLSApprox::OpMWLSRhoAtGaussUsingPrecalulatedCoeffs::calculate2ndDerivative
const bool calculate2ndDerivative
Definition: MWLS.hpp:307
FractureMechanics
Definition: AnalyticalFun.hpp:15
FractureMechanics::MWLSApprox::distLeafs
std::vector< std::pair< double, EntityHandle > > distLeafs
Definition: MWLS.hpp:826
FractureMechanics::MWLSApprox::OpMWLSMaterialStressLhs_Dx::doWork
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
Definition: MWLS.cpp:2065
MoFEM::PetscOptionsGetBool
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
Definition: DeprecatedPetsc.hpp:182
Generated by Doxygen 1.8.17 and hosted at University of Glasgow