v0.14.0
CrackFrontElement.cpp
Go to the documentation of this file.
1 /** \file CrackFrontElement.cpp
2  */
3 
4 /* This file is part of MoFEM.
5  * MoFEM is free software: you can redistribute it and/or modify it under
6  * the terms of the GNU Lesser General Public License as published by the
7  * Free Software Foundation, either version 3 of the License, or (at your
8  * option) any later version.
9  *
10  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
11  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
12  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
13  * License for more details.
14  *
15  * You should have received a copy of the GNU Lesser General Public
16  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
17 
18 #include <moab/SpatialLocator.hpp>
19 #include <moab/ElemEvaluator.hpp>
20 
21 #include <MoFEM.hpp>
22 using namespace MoFEM;
23 
24 #include <cholesky.hpp>
25 
26 #include <BasicFiniteElements.hpp>
27 #include <Mortar.hpp>
28 #include <Hooke.hpp>
29 
30 
31 #include <AnalyticalFun.hpp>
32 
33 
34 // TODO implementation in ComplexConstArea is obsolete, not follows naming
35 // convention and code is unnecessarily complicated and in places difficult to
36 // follow.
37 #include <ComplexConstArea.hpp>
38 #include <ConstantArea.hpp>
39 
40 extern "C" {
41 #include <phg-quadrule/quad.h>
42 }
43 
44 #include <NeoHookean.hpp>
45 #include <Hooke.hpp>
46 
47 #include <CrackFrontElement.hpp>
48 
49 namespace po = boost::program_options;
50 
51 namespace FractureMechanics {
52 
53 static std::map<long int, MatrixDouble> mapRefCoords;
54 
55 template <>
56 template <>
57 int CrackFrontSingularBase<NonlinearElasticElement::MyVolumeFE,
58  VolumeElementForcesAndSourcesCore>::
59  getRuleImpl(TwoType<NonlinearElasticElement::MyVolumeFE,
60  VolumeElementForcesAndSourcesCore>,
61  int order) {
62  return -1;
63 }
64 
65 template <>
66 template <>
67 MoFEMErrorCode CrackFrontSingularBase<NonlinearElasticElement::MyVolumeFE,
68  VolumeElementForcesAndSourcesCore>::
69  setGaussPtsImpl(TwoType<NonlinearElasticElement::MyVolumeFE,
70  VolumeElementForcesAndSourcesCore>,
71  int order) {
72 
73  MoFEMFunctionBeginHot;
74 
75  auto set_gauss_pts = [this]() {
76  MoFEMFunctionBegin;
77  gaussPts.swap(refGaussCoords);
78  const size_t nb_gauss_pts = gaussPts.size2();
79  dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE).resize(nb_gauss_pts, 4);
80  double *shape_ptr =
81  &*dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE).data().begin();
82  CHKERR ShapeMBTET(shape_ptr, &gaussPts(0, 0), &gaussPts(1, 0),
83  &gaussPts(2, 0), nb_gauss_pts);
84  MoFEMFunctionReturn(0);
85  };
86 
87  std::bitset<4> singular_nodes_ids;
88  if (refineIntegration == PETSC_TRUE && singularElement) {
89 
90  for (int nn = 0; nn != 4; nn++) {
91  if (singularNodes[nn])
92  singular_nodes_ids.set(nn);
93  }
94 
95  auto it_map_ref_coords = mapRefCoords.find(singular_nodes_ids.to_ulong());
96  if (it_map_ref_coords != mapRefCoords.end()) {
97 
98  refGaussCoords = it_map_ref_coords->second;
99  CHKERR set_gauss_pts();
100 
101  MoFEMFunctionReturnHot(0);
102  }
103  }
104 
105  order = order ? order : 1;
106  int rule = (singularElement) ? 2 * (std::max(order, 3) - 1) : 2 * (order - 1);
107 
108  if (rule < QUAD_3D_TABLE_SIZE) {
109  if (QUAD_3D_TABLE[rule]->dim != 3) {
110  SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "wrong dimension");
111  }
112  if (QUAD_3D_TABLE[rule]->order < rule) {
113  SETERRQ2(mField.get_comm(), MOFEM_DATA_INCONSISTENCY,
114  "wrong order %d != %d", QUAD_3D_TABLE[rule]->order, rule);
115  }
116  const size_t nb_gauss_pts = QUAD_3D_TABLE[rule]->npoints;
117  gaussPts.resize(4, nb_gauss_pts, false);
118  cblas_dcopy(nb_gauss_pts, &QUAD_3D_TABLE[rule]->points[1], 4, &gaussPts(0, 0),
119  1);
120  cblas_dcopy(nb_gauss_pts, &QUAD_3D_TABLE[rule]->points[2], 4, &gaussPts(1, 0),
121  1);
122  cblas_dcopy(nb_gauss_pts, &QUAD_3D_TABLE[rule]->points[3], 4, &gaussPts(2, 0),
123  1);
124  cblas_dcopy(nb_gauss_pts, QUAD_3D_TABLE[rule]->weights, 1, &gaussPts(3, 0),
125  1);
126  dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE).resize(nb_gauss_pts, 4,
127  false);
128  double *shape_ptr =
129  &*dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE).data().begin();
130  cblas_dcopy(4 * nb_gauss_pts, QUAD_3D_TABLE[rule]->points, 1, shape_ptr, 1);
131  } else {
132  SETERRQ2(mField.get_comm(), MOFEM_DATA_INCONSISTENCY,
133  "rule > quadrature order %d < %d", rule, QUAD_3D_TABLE_SIZE);
134  }
135 
136  // Refine quadrature by splitting edges adjacent to crack front. On sub
137  // triangle place quadrature points.
138  if (!singularElement)
139  MoFEMFunctionReturnHot(0);
140 
141  // refine mesh on refine elements
142  if (refineIntegration == PETSC_TRUE) {
143 
144  const int max_level = refinementLevels;
145  EntityHandle tet;
146 
147  moab::Core moab_ref;
148  double base_coords[] = {0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1};
149  EntityHandle nodes[4];
150  for (int nn = 0; nn != 4; nn++)
151  CHKERR moab_ref.create_vertex(&base_coords[3 * nn], nodes[nn]);
152  CHKERR moab_ref.create_element(MBTET, nodes, 4, tet);
153  MoFEM::CoreTmp<-1> mofem_ref_core(moab_ref, PETSC_COMM_SELF, -2);
154  MoFEM::Interface &m_field_ref = mofem_ref_core;
155  CHKERR m_field_ref.getInterface<BitRefManager>()->setBitRefLevelByDim(
156  0, 3, BitRefLevel().set(0));
157 
158  Range singular_nodes;
159  for (int nn = 0; nn != 4; nn++) {
160  if (singularNodes[nn]) {
161  EntityHandle ent;
162  CHKERR moab_ref.side_element(tet, 0, nn, ent);
163  singular_nodes.insert(ent);
164  }
165  }
166 
167  EntityHandle meshset;
168  if (singular_nodes.size() > 1) {
169  CHKERR moab_ref.create_meshset(MESHSET_SET, meshset);
170  for (int ee = 0; ee != 6; ee++) {
171  if (singularEdges[ee]) {
172  EntityHandle ent;
173  CHKERR moab_ref.side_element(tet, 1, ee, ent);
174  CHKERR moab_ref.add_entities(meshset, &ent, 1);
175  }
176  }
177  }
178 
179  MeshRefinement *m_ref;
180  CHKERR m_field_ref.getInterface(m_ref);
181  for (int ll = 0; ll != max_level; ll++) {
182  // PetscPrintf(T::mField.get_comm(),"Refine Level %d\n",ll);
183  Range edges;
184  CHKERR m_field_ref.getInterface<BitRefManager>()
185  ->getEntitiesByTypeAndRefLevel(BitRefLevel().set(ll),
186  BitRefLevel().set(), MBEDGE, edges);
187  Range ref_edges;
188  CHKERR moab_ref.get_adjacencies(singular_nodes, 1, true, ref_edges,
189  moab::Interface::UNION);
190  ref_edges = intersect(ref_edges, edges);
191  if (singular_nodes.size() > 1) {
192  Range ents;
193  CHKERR moab_ref.get_entities_by_type(meshset, MBEDGE, ents, true);
194  ref_edges = intersect(ref_edges, ents);
195  }
196  Range tets;
197  CHKERR m_field_ref.getInterface<BitRefManager>()
198  ->getEntitiesByTypeAndRefLevel(BitRefLevel().set(ll),
199  BitRefLevel().set(), MBTET, tets);
200  // refine mesh
201  CHKERR m_ref->addVerticesInTheMiddleOfEdges(
202  ref_edges, BitRefLevel().set(ll + 1));
203  CHKERR m_ref->refineTets(tets, BitRefLevel().set(ll + 1));
204  if (singular_nodes.size() > 1) {
205  CHKERR m_field_ref.getInterface<BitRefManager>()
206  ->updateMeshsetByEntitiesChildren(
207  meshset, BitRefLevel().set(ll + 1), meshset, MBEDGE, true);
208  }
209  }
210 
211  Range tets;
212  CHKERR m_field_ref.getInterface<BitRefManager>()
213  ->getEntitiesByTypeAndRefLevel(BitRefLevel().set(max_level),
214  BitRefLevel().set(), MBTET, tets);
215  refCoords.resize(tets.size(), 12, false);
216  int tt = 0;
217  for (Range::iterator tit = tets.begin(); tit != tets.end(); tit++, tt++) {
218  int num_nodes;
219  const EntityHandle *conn;
220  CHKERR moab_ref.get_connectivity(*tit, conn, num_nodes, false);
221  CHKERR moab_ref.get_coords(conn, num_nodes, &refCoords(tt, 0));
222  }
223 
224  if (0) {
225  EntityHandle out_meshset;
226  CHKERR moab_ref.create_meshset(MESHSET_SET, out_meshset);
227  CHKERR moab_ref.add_entities(out_meshset, tets);
228  std::string file_name =
229  "ref_tet_" + boost::lexical_cast<std::string>(nInTheLoop) + ".vtk";
230  CHKERR moab_ref.write_file(file_name.c_str(), "VTK", "", &out_meshset, 1);
231  CHKERR moab_ref.delete_entities(&out_meshset, 1);
232  }
233 
234  // Set integration points
235  double diff_n[12];
236  CHKERR ShapeDiffMBTET(diff_n);
237  const size_t nb_gauss_pts = gaussPts.size2();
238  refGaussCoords.resize(4, nb_gauss_pts * refCoords.size1());
239  MatrixDouble &shape_n = dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE);
240  int gg = 0;
241  for (size_t tt = 0; tt != refCoords.size1(); tt++) {
242  double *tet_coords = &refCoords(tt, 0);
243  double det = ShapeVolumeMBTET(diff_n, tet_coords);
244  det *= 6;
245  for (size_t ggg = 0; ggg != nb_gauss_pts; ++ggg, ++gg) {
246  for (int dd = 0; dd != 3; dd++) {
247  refGaussCoords(dd, gg) = shape_n(ggg, 0) * tet_coords[3 * 0 + dd] +
248  shape_n(ggg, 1) * tet_coords[3 * 1 + dd] +
249  shape_n(ggg, 2) * tet_coords[3 * 2 + dd] +
250  shape_n(ggg, 3) * tet_coords[3 * 3 + dd];
251  }
252  refGaussCoords(3, gg) = gaussPts(3, ggg) * det;
253  }
254  }
255 
256  mapRefCoords[singular_nodes_ids.to_ulong()] = refGaussCoords;
257  CHKERR set_gauss_pts();
258  }
259 
260  MoFEMFunctionReturnHot(0);
261 }
262 
263 MoFEMErrorCode OpGetCrackFrontDataGradientAtGaussPts::doWork(
264  int side, EntityType type, DataForcesAndSourcesCore::EntData &data) {
265  MoFEMFunctionBegin;
266 
267  if (!singularElement) {
268  return OpCalculateVectorFieldGradient<3, 3>::doWork(side, type, data);
269  }
270 
271  const int nb_dofs = data.getFieldData().size();
272  if (!nb_dofs && type == this->zeroType) {
273  this->dataPtr->resize(9, 0, false);
274  MoFEMFunctionReturnHot(0);
275  }
276  if (!nb_dofs) {
277  MoFEMFunctionReturnHot(0);
278  }
279  const int nb_gauss_pts = data.getN().size1();
280  const int nb_base_functions = data.getN().size2();
281  ublas::matrix<double, ublas::row_major, DoubleAllocator> &mat =
282  *this->dataPtr;
283  if (type == this->zeroType) {
284  mat.resize(9, nb_gauss_pts, false);
285  mat.clear();
286  }
287  auto diff_base_function = data.getFTensor1DiffN<3>();
288  auto gradients_at_gauss_pts = getFTensor2FromMat<3, 3>(mat);
289  FTensor::Index<'I', 3> I;
290  FTensor::Index<'J', 3> J;
291  FTensor::Index<'K', 3> K;
292  int size = nb_dofs / 3;
293  if (nb_dofs % 3) {
294  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Data inconsistency");
295  }
296  double *inv_jac_ptr = &*invSJac.data().begin();
297  FTensor::Tensor2<double *, 3, 3> t_inv_singular_jacobian(
298  inv_jac_ptr, &inv_jac_ptr[1], &inv_jac_ptr[2], &inv_jac_ptr[3],
299  &inv_jac_ptr[4], &inv_jac_ptr[5], &inv_jac_ptr[6], &inv_jac_ptr[7],
300  &inv_jac_ptr[8], 9);
301 
302  for (int gg = 0; gg < nb_gauss_pts; ++gg) {
303  auto field_data = data.getFTensor1FieldData<3>();
304  int bb = 0;
305  for (; bb < size; ++bb) {
306  gradients_at_gauss_pts(I, J) +=
307  field_data(I) *
308  (t_inv_singular_jacobian(K, J) * diff_base_function(K));
309  ++field_data;
310  ++diff_base_function;
311  }
312  // Number of dofs can be smaller than number of Tensor_Dim0 x base functions
313  for (; bb != nb_base_functions; ++bb)
314  ++diff_base_function;
315  ++gradients_at_gauss_pts;
316  ++t_inv_singular_jacobian;
317  }
318 
319  MoFEMFunctionReturn(0);
320 }
321 
322 MoFEMErrorCode
323 OpGetCrackFrontDataAtGaussPts::doWork(int side, EntityType type,
324  DataForcesAndSourcesCore::EntData &data) {
325  MoFEMFunctionBegin;
326 
327  if (!singularElement) {
328  return NonlinearElasticElement::OpGetDataAtGaussPts::doWork(side, type,
329  data);
330  }
331 
332  const int nb_dofs = data.getFieldData().size();
333  const int nb_base_functions = data.getN().size2();
334  if (nb_dofs == 0) {
335  MoFEMFunctionReturnHot(0);
336  }
337  const int nb_gauss_pts = data.getN().size1();
338  const int rank = data.getFieldDofs()[0]->getNbOfCoeffs();
339 
340  // initialize
341  if (type == zeroAtType) {
342  valuesAtGaussPts.resize(nb_gauss_pts);
343  gradientAtGaussPts.resize(nb_gauss_pts);
344  for (int gg = 0; gg != nb_gauss_pts; gg++) {
345  valuesAtGaussPts[gg].resize(rank, false);
346  gradientAtGaussPts[gg].resize(rank, 3, false);
347  }
348  for (int gg = 0; gg != nb_gauss_pts; gg++) {
349  valuesAtGaussPts[gg].clear();
350  gradientAtGaussPts[gg].clear();
351  }
352  }
353 
354  auto base_function = data.getFTensor0N();
355  auto diff_base_functions = data.getFTensor1DiffN<3>();
356  FTensor::Index<'i', 3> i;
357  FTensor::Index<'j', 3> j;
358  FTensor::Index<'k', 3> k;
359 
360  double *inv_jac_ptr = &*invSJac.data().begin();
361  FTensor::Tensor2<double *, 3, 3> t_inv_singular_jacobian(
362  inv_jac_ptr, &inv_jac_ptr[1], &inv_jac_ptr[2], &inv_jac_ptr[3],
363  &inv_jac_ptr[4], &inv_jac_ptr[5], &inv_jac_ptr[6], &inv_jac_ptr[7],
364  &inv_jac_ptr[8], 9);
365 
366  if (rank == 3) {
367 
368  for (int gg = 0; gg != nb_gauss_pts; gg++) {
369  auto field_data = data.getFTensor1FieldData<3>();
370  FTensor::Tensor1<double *, 3> values(&valuesAtGaussPts[gg][0],
371  &valuesAtGaussPts[gg][1],
372  &valuesAtGaussPts[gg][2]);
373  FTensor::Tensor2<double *, 3, 3> gradient(
374  &gradientAtGaussPts[gg](0, 0), &gradientAtGaussPts[gg](0, 1),
375  &gradientAtGaussPts[gg](0, 2), &gradientAtGaussPts[gg](1, 0),
376  &gradientAtGaussPts[gg](1, 1), &gradientAtGaussPts[gg](1, 2),
377  &gradientAtGaussPts[gg](2, 0), &gradientAtGaussPts[gg](2, 1),
378  &gradientAtGaussPts[gg](2, 2));
379  int bb = 0;
380  for (; bb != nb_dofs / 3; bb++) {
381  values(i) += base_function * field_data(i);
382  gradient(i, j) += field_data(i) * (t_inv_singular_jacobian(k, j) *
383  diff_base_functions(k));
384  ++diff_base_functions;
385  ++base_function;
386  ++field_data;
387  }
388  for (; bb != nb_base_functions; bb++) {
389  ++diff_base_functions;
390  ++base_function;
391  }
392  ++t_inv_singular_jacobian;
393  }
394 
395  } else {
396 
397  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
398  "Not implemented for rank = %d", rank);
399  }
400 
401  MoFEMFunctionReturn(0);
402 }
403 
404 MoFEMErrorCode
405 OpPostProcDisplacements::doWork(int side, EntityType type,
406  DataForcesAndSourcesCore::EntData &data) {
407  MoFEMFunctionBegin;
408  if (!singularElement || type != MBVERTEX) {
409  MoFEMFunctionReturnHot(0);
410  }
411  double def_VAL[3] = {0, 0, 0};
412  Tag th_singular_disp;
413  CHKERR postProcMesh.tag_get_handle("SINGULAR_DISP", 3, MB_TYPE_DOUBLE,
414  th_singular_disp,
415  MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
416  FTensor::Index<'i', 3> i;
417  FTensor::Index<'j', 3> j;
418  FTensor::Index<'k', 3> k;
419  auto t_H = getFTensor2FromMat<3, 3>(*H);
420  FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_singular_disp(
421  &singularDisp(0, 0), &singularDisp(0, 1), &singularDisp(0, 2));
422  FTensor::Tensor1<double, 3> t_singular_current_disp;
423  for (unsigned int gg = 0; gg != singularDisp.size1(); gg++) {
424  t_singular_current_disp(i) = t_H(i, j) * t_singular_disp(j);
425  CHKERR postProcMesh.tag_set_data(th_singular_disp, &mapGaussPts[gg], 1,
426  &t_singular_current_disp(0));
427  ++t_H;
428  ++t_singular_disp;
429  }
430  MoFEMFunctionReturn(0);
431 }
432 
433 MoFEMErrorCode OpTransfromSingularBaseFunctions::doWork(
434  int side, EntityType type, DataForcesAndSourcesCore::EntData &data) {
435  FTensor::Index<'i', 3> i;
436  FTensor::Index<'j', 3> j;
437  MoFEMFunctionBegin;
438 
439  if (!singularElement)
440  MoFEMFunctionReturnHot(0);
441 
442  if (invSJac.size2() != 9) {
443  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
444  "It looks that ho inverse of Jacobian is not calculated %d != 9",
445  invSJac.size2());
446  }
447  double *t_inv_jac_ptr = &*invSJac.data().begin();
448  FTensor::Tensor2<double *, 3, 3> t_inv_jac(
449  t_inv_jac_ptr, &t_inv_jac_ptr[1], &t_inv_jac_ptr[2], &t_inv_jac_ptr[3],
450  &t_inv_jac_ptr[4], &t_inv_jac_ptr[5], &t_inv_jac_ptr[6],
451  &t_inv_jac_ptr[7], &t_inv_jac_ptr[8], 9);
452 
453  switch (type) {
454  case MBVERTEX:
455  if (applyDet) {
456  for (unsigned int gg = 0; gg != getGaussPts().size2(); ++gg) {
457  getGaussPts()(3, gg) *= detS[gg];
458  }
459  }
460  case MBEDGE:
461  case MBTRI:
462  case MBTET: {
463  unsigned int nb_gauss_pts = data.getN(bAse).size1();
464  if (nb_gauss_pts == 0)
465  MoFEMFunctionReturnHot(0);
466  if (invSJac.size1() != nb_gauss_pts) {
467  SETERRQ2(
468  PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
469  "It looks that ho inverse of Jacobian is not calculated %d != %d",
470  invSJac.size1(), nb_gauss_pts);
471  }
472  unsigned int nb_base_functions = data.getN(bAse).size2();
473  if (nb_base_functions == 0)
474  MoFEMFunctionReturnHot(0);
475  diffNinvJac.resize(nb_gauss_pts, 3 * nb_base_functions, false);
476  double *t_inv_n_ptr = &*diffNinvJac.data().begin();
477  FTensor::Tensor1<double *, 3> t_inv_diff_n(t_inv_n_ptr, &t_inv_n_ptr[1],
478  &t_inv_n_ptr[2], 3);
479  // Apply transformation to each bAse function ob given element
480  FTensor::Tensor1<double *, 3> t_diff_n = FTensor::Tensor1<double *, 3>(
481  &data.getDiffN(bAse)(0, 0), &data.getDiffN(bAse)(0, 1),
482  &data.getDiffN(bAse)(0, 2), 3);
483  for (unsigned int gg = 0; gg != nb_gauss_pts; ++gg) {
484  for (unsigned int bb = 0; bb != nb_base_functions; ++bb) {
485  t_inv_diff_n(i) = t_diff_n(j) * t_inv_jac(j, i);
486  ++t_diff_n;
487  ++t_inv_diff_n;
488  }
489  ++t_inv_jac;
490  }
491  if (type == MBVERTEX) {
492  data.getDiffN(bAse).resize(diffNinvJac.size1(), diffNinvJac.size2(),
493  false);
494  }
495  data.getDiffN(bAse).data().swap(diffNinvJac.data());
496  } break;
497  default:
498  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "not implemented");
499  }
500 
501  MoFEMFunctionReturn(0);
502 }
503 
504 OpRhsBoneExplicitDerivariveWithHooke::OpRhsBoneExplicitDerivariveWithHooke(
505  HookeElement::DataAtIntegrationPts &common_data,
506  boost::shared_ptr<VectorDouble> rho_at_gauss_pts,
507  boost::shared_ptr<MatrixDouble> rho_grad_at_gauss_pts, Range tets_in_block,
508  const double rho_0, const double bone_n, Range *forces_only_on_entities_row)
509  : MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator(
510  "MESH_NODE_POSITIONS", UserDataOperator::OPROW),
511  commonData(common_data), rhoAtGaussPts(rho_at_gauss_pts),
512  rhoGradAtGaussPts(rho_grad_at_gauss_pts), tetsInBlock(tets_in_block),
513  rHo0(rho_0), boneN(bone_n) {
514  if (forces_only_on_entities_row)
515  forcesOnlyOnEntitiesRow = *forces_only_on_entities_row;
516 }
517 
518 MoFEMErrorCode OpRhsBoneExplicitDerivariveWithHooke::doWork(
519  int row_side, EntityType row_type,
520  DataForcesAndSourcesCore::EntData &row_data) {
521  FTensor::Index<'i', 3> i;
522  MoFEMFunctionBeginHot;
523  const int nb_dofs = row_data.getIndices().size();
524  if (!nb_dofs)
525  MoFEMFunctionReturnHot(0);
526  if (tetsInBlock.find(getFEEntityHandle()) == tetsInBlock.end()) {
527  MoFEMFunctionReturnHot(0);
528  }
529  // If this is not part of blockset
530  const int nb_gauss_pts = row_data.getN().size1();
531  nF.resize(nb_dofs, false);
532  nF.clear();
533  auto t_rho_gradient = getFTensor1FromMat<3>(*rhoGradAtGaussPts);
534  auto rho = getFTensor0FromVec(*rhoAtGaussPts);
535  auto energy = getFTensor0FromVec(*commonData.energyVec);
536  auto det_H = getFTensor0FromVec(*commonData.detHVec);
537 
538  for (int gg = 0; gg != nb_gauss_pts; gg++) {
539  const double w = getVolume() * getGaussPts()(3, gg) * det_H;
540 
541  const double denergy_drho =
542  energy * (boneN / rho); // derivative of energy w.r.t density
543  auto t_base = row_data.getFTensor0N(gg, 0);
544  FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_nf(&nF[0], &nF[1],
545  &nF[2]);
546  for (int bb = 0; bb != nb_dofs / 3; bb++) {
547  t_nf(i) += t_base * (w * denergy_drho * t_rho_gradient(i));
548  ++t_nf;
549  ++t_base;
550  }
551  ++t_rho_gradient;
552  ++rho;
553  ++energy;
554  ++det_H;
555  }
556  int *indices_ptr = &row_data.getIndices()[0];
557  if (!forcesOnlyOnEntitiesRow.empty()) {
558  iNdices.resize(nb_dofs, false);
559  noalias(iNdices) = row_data.getIndices();
560  indices_ptr = &iNdices[0];
561  auto &dofs = row_data.getFieldDofs();
562  auto dit = dofs.begin();
563  for (int ii = 0; dit != dofs.end(); dit++, ii++) {
564  if (auto dof = (*dit)) {
565  if (forcesOnlyOnEntitiesRow.find(dof->getEnt()) ==
566  forcesOnlyOnEntitiesRow.end()) {
567  iNdices[ii] = -1;
568  }
569  }
570  }
571  }
572 
573  CHKERR VecSetOption(getFEMethod()->snes_f, VEC_IGNORE_NEGATIVE_INDICES,
574  PETSC_TRUE);
575  CHKERR VecSetValues(getFEMethod()->snes_f, nb_dofs, indices_ptr, &nF[0],
576  ADD_VALUES);
577  MoFEMFunctionReturnHot(0);
578 }
579 
580 OpLhsBoneExplicitDerivariveWithHooke_dX::
581  OpLhsBoneExplicitDerivariveWithHooke_dX(
582  HookeElement::DataAtIntegrationPts &common_data,
583  boost::shared_ptr<VectorDouble> rho_at_gauss_pts,
584  boost::shared_ptr<MatrixDouble> diff_rho_at_gauss_pts,
585  boost::shared_ptr<MatrixDouble> diff_diff_rho_at_gauss_pts,
586  MatrixDouble &singular_displ, Range tets_in_block, const double rho_0,
587  const double bone_n, Range *forces_only_on_entities_row)
588  : MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator(
589  "MESH_NODE_POSITIONS", "MESH_NODE_POSITIONS",
590  UserDataOperator::OPROWCOL, false), // FIXME:
591  commonData(common_data), rhoAtGaussPts(rho_at_gauss_pts),
592  diffRhoAtGaussPts(diff_rho_at_gauss_pts),
593  diffDiffRhoAtGaussPts(diff_diff_rho_at_gauss_pts),
594  singularDispl(singular_displ), tetsInBlock(tets_in_block), rHo0(rho_0),
595  boneN(bone_n) {
596  if (forces_only_on_entities_row)
597  forcesOnlyOnEntitiesRow = *forces_only_on_entities_row;
598 }
599 
600 MoFEMErrorCode OpLhsBoneExplicitDerivariveWithHooke_dX::doWork(
601  int row_side, int col_side, EntityType row_type, EntityType col_type,
602  DataForcesAndSourcesCore::EntData &row_data,
603  DataForcesAndSourcesCore::EntData &col_data) {
604  MoFEMFunctionBegin;
605  if (tetsInBlock.find(getFEEntityHandle()) == tetsInBlock.end()) {
606  MoFEMFunctionReturnHot(0);
607  }
608 
609  const int row_nb_dofs = row_data.getIndices().size();
610  if (!row_nb_dofs)
611  MoFEMFunctionReturnHot(0);
612 
613  const int col_nb_dofs = col_data.getIndices().size();
614  if (!col_nb_dofs)
615  MoFEMFunctionReturnHot(0);
616 
617  nA.resize(row_nb_dofs, col_nb_dofs, false);
618  nA.clear();
619 
620  const int nb_gauss_pts = row_data.getN().size1();
621  auto t_invH = getFTensor2FromMat<3, 3>(*commonData.invHMat);
622  auto t_h = getFTensor2FromMat<3, 3>(*commonData.hMat);
623  auto t_H = getFTensor2FromMat<3, 3>(*commonData.HMat);
624  auto rho = getFTensor0FromVec(*rhoAtGaussPts);
625  auto t_diff_rho = getFTensor1FromMat<3>(*diffRhoAtGaussPts);
626  auto t_diff_diff_rho = getFTensor2FromMat<3, 3>(*diffDiffRhoAtGaussPts);
627 
628  auto t_cauchy_stress =
629  getFTensor2SymmetricFromMat<3>(*commonData.cauchyStressMat);
630 
631  auto energy = getFTensor0FromVec(*commonData.energyVec);
632  auto det_H = getFTensor0FromVec(*commonData.detHVec);
633 
634  auto t_initial_singular_displacement = getFTensor1FromMat<3>(singularDispl);
635  auto t_row_diff_n = row_data.getFTensor1DiffN<3>();
636 
637  const int nb_row_base = row_data.getN().size2();
638 
639  FTensor::Index<'i', 3> i;
640  FTensor::Index<'j', 3> j;
641  FTensor::Index<'k', 3> k;
642  FTensor::Index<'l', 3> l;
643  FTensor::Index<'m', 3> m;
644  FTensor::Index<'n', 3> n;
645 
646  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
647  const double val = getVolume() * getGaussPts()(3, gg) * det_H;
648 
649  const double drho_coeff = (boneN / rho);
650  const double denergy_drho = energy * drho_coeff;
651  const double denergy_drho_drho =
652  energy * (boneN - 1.) * boneN / (rho * rho);
653  FTensor::Tensor4<double, 3, 3, 3, 3> t_F_dX;
654  t_F_dX(i, j, k, l) = -(t_h(i, m) * t_invH(m, k)) * t_H(l, j);
655 
656  FTensor::Tensor2<double, 3, 3> t_energy_dX;
657  t_energy_dX(k, l) = t_F_dX(i, j, k, l) * t_cauchy_stress(i, j);
658 
659  FTensor::Tensor1<double, 3> t_b;
660  FTensor::Tensor1<double, 3> t_rho_dX_pulled;
661  FTensor::Tensor1<double, 3> t_psi_dX_rho_dX;
662  FTensor::Tensor1<double, 3> t_psi_dX_rho_dX_pulled;
663  FTensor::Tensor2<double, 3, 3> t_rho_dX_dX_pulled;
664  FTensor::Tensor2<double, 3, 3> t_energy_dX_pulled;
665  t_rho_dX_pulled(i) = t_diff_rho(j) * t_invH(j, i);
666  t_psi_dX_rho_dX(i) =
667  t_energy_dX(i, j) * t_rho_dX_pulled(j) * val * drho_coeff;
668  t_psi_dX_rho_dX_pulled(j) = t_psi_dX_rho_dX(i) * t_invH(i, j);
669 
670  auto t_row_base = row_data.getFTensor0N(gg, 0);
671  int rr = 0;
672  for (; rr != row_nb_dofs / 3; rr++) {
673  FTensor::Tensor2<double *, 3, 3> t_mat(
674  &nA(3 * rr + 0, 0), &nA(3 * rr + 0, 1), &nA(3 * rr + 0, 2),
675  &nA(3 * rr + 1, 0), &nA(3 * rr + 1, 1), &nA(3 * rr + 1, 2),
676  &nA(3 * rr + 2, 0), &nA(3 * rr + 2, 1), &nA(3 * rr + 2, 2), 3);
677 
678  FTensor::Tensor2<double, 3, 3> t_a;
679  t_a(i, j) = val * denergy_drho_drho * t_diff_rho(i) * t_diff_rho(j);
680  t_a(i, j) += val * denergy_drho * t_diff_diff_rho(i, j);
681  t_b(i) = t_row_base * val * denergy_drho * t_diff_rho(i);
682 
683  auto t_col_base = col_data.getFTensor0N(gg, 0);
684  auto t_col_diff_base = col_data.getFTensor1DiffN<3>(gg, 0);
685 
686  for (int cc = 0; cc != col_nb_dofs / 3; cc++) {
687 
688  t_mat(i, j) += t_row_base * t_a(i, j) * t_col_base;
689  t_mat(i, j) += t_row_base * t_a(i, j) * t_col_diff_base(k) *
690  t_initial_singular_displacement(k);
691 
692  t_mat(i, j) += t_b(i) * (t_invH(l, j) * t_col_diff_base(l));
693  t_mat(i, j) +=
694  t_row_base * t_psi_dX_rho_dX_pulled(i) * t_col_diff_base(j);
695  ++t_mat;
696  ++t_col_base;
697  ++t_col_diff_base;
698  }
699  ++t_row_base;
700  ++t_row_diff_n;
701  }
702  for (; rr != nb_row_base; ++rr) {
703  ++t_row_base;
704  ++t_row_diff_n;
705  }
706 
707  ++t_invH;
708  ++t_H;
709  ++rho;
710  ++energy;
711  ++det_H;
712  ++t_diff_diff_rho;
713  ++t_diff_rho;
714  ++t_h;
715  // ++t_D;
716  // ++t_strain;
717  ++t_cauchy_stress;
718  ++t_initial_singular_displacement;
719  }
720 
721  // Assembly
722  int *indices_ptr = &row_data.getIndices()[0];
723  if (!forcesOnlyOnEntitiesRow.empty()) {
724  iNdices.resize(row_nb_dofs);
725  noalias(iNdices) = row_data.getIndices();
726  indices_ptr = &iNdices[0];
727  auto &dofs = row_data.getFieldDofs();
728  auto dit = dofs.begin();
729  for (int ii = 0; dit != dofs.end(); dit++, ii++) {
730  if (auto dof = (*dit)) {
731  if (forcesOnlyOnEntitiesRow.find(dof->getEnt()) ==
732  forcesOnlyOnEntitiesRow.end()) {
733  iNdices[ii] = -1;
734  }
735  }
736  }
737  }
738 
739  CHKERR MatSetValues(getFEMethod()->snes_B, row_nb_dofs, indices_ptr,
740  col_nb_dofs, &*col_data.getIndices().begin(),
741  &*nA.data().begin(), ADD_VALUES);
742  MoFEMFunctionReturn(0);
743 }
744 
745 OpLhsBoneExplicitDerivariveWithHooke_dx::
746  OpLhsBoneExplicitDerivariveWithHooke_dx(
747  HookeElement::DataAtIntegrationPts &common_data,
748  boost::shared_ptr<VectorDouble> rho_at_gauss_pts,
749  boost::shared_ptr<MatrixDouble> diff_rho_at_gauss_pts,
750  boost::shared_ptr<MatrixDouble> diff_diff_rho_at_gauss_pts,
751  MatrixDouble &singular_displ, Range tets_in_block, const double rho_0,
752  const double bone_n, Range *forces_only_on_entities_row)
753  : MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator(
754  "MESH_NODE_POSITIONS", "SPATIAL_POSITION", UserDataOperator::OPROWCOL,
755  false),
756  commonData(common_data), rhoAtGaussPts(rho_at_gauss_pts),
757  diffRhoAtGaussPts(diff_rho_at_gauss_pts),
758  diffDiffRhoAtGaussPts(diff_diff_rho_at_gauss_pts),
759  singularDispl(singular_displ), tetsInBlock(tets_in_block), rHo0(rho_0),
760  boneN(bone_n) {
761  if (forces_only_on_entities_row)
762  forcesOnlyOnEntitiesRow = *forces_only_on_entities_row;
763 }
764 
765 MoFEMErrorCode OpLhsBoneExplicitDerivariveWithHooke_dx::doWork(
766  int row_side, int col_side, EntityType row_type, EntityType col_type,
767  DataForcesAndSourcesCore::EntData &row_data,
768  DataForcesAndSourcesCore::EntData &col_data) {
769  MoFEMFunctionBegin;
770  if (tetsInBlock.find(getFEEntityHandle()) == tetsInBlock.end()) {
771  MoFEMFunctionReturnHot(0);
772  }
773  const int row_nb_dofs = row_data.getIndices().size();
774  if (!row_nb_dofs) {
775  MoFEMFunctionReturnHot(0);
776  }
777  const int col_nb_dofs = col_data.getIndices().size();
778  if (!col_nb_dofs) {
779  MoFEMFunctionReturnHot(0);
780  }
781  nA.resize(row_nb_dofs, col_nb_dofs, false);
782  nA.clear();
783 
784  const int nb_row_base = row_data.getN().size2();
785  const int nb_gauss_pts = row_data.getN().size1();
786 
787  auto t_initial_singular_displacement = getFTensor1FromMat<3>(singularDispl);
788  auto t_row_diff_n = row_data.getFTensor1DiffN<3>();
789 
790  auto t_H = getFTensor2FromMat<3, 3>(*(commonData.HMat));
791  auto t_h = getFTensor2FromMat<3, 3>(*(commonData.hMat));
792  auto det_H = getFTensor0FromVec(*commonData.detHVec);
793  auto t_invH = getFTensor2FromMat<3, 3>(*commonData.invHMat);
794  auto energy = getFTensor0FromVec(*commonData.energyVec);
795 
796  auto t_cauchy_stress =
797  getFTensor2SymmetricFromMat<3>(*(commonData.cauchyStressMat));
798  auto rho = getFTensor0FromVec(*rhoAtGaussPts);
799 
800  auto t_diff_rho = getFTensor1FromMat<3>(*diffRhoAtGaussPts);
801 
802  FTensor::Index<'i', 3> i;
803  FTensor::Index<'j', 3> j;
804  FTensor::Index<'k', 3> k;
805  FTensor::Index<'l', 3> l;
806  FTensor::Index<'m', 3> m;
807  FTensor::Tensor1<double, 3> t_a;
808 
809  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
810  double val = getVolume() * getGaussPts()(3, gg) * det_H;
811 
812  auto t_row_base = row_data.getFTensor0N(gg, 0);
813  const double drho_coef = (boneN / rho);
814 
815  int rr = 0;
816  for (; rr != row_nb_dofs / 3; ++rr) {
817 
818  auto t_col_n = col_data.getFTensor0N(gg, 0);
819  auto t_col_diff_base = col_data.getFTensor1DiffN<3>(gg, 0);
820 
821  FTensor::Tensor2<double *, 3, 3> t_mat(
822  &nA(3 * rr + 0, 0), &nA(3 * rr + 0, 1), &nA(3 * rr + 0, 2),
823  &nA(3 * rr + 1, 0), &nA(3 * rr + 1, 1), &nA(3 * rr + 1, 2),
824  &nA(3 * rr + 2, 0), &nA(3 * rr + 2, 1), &nA(3 * rr + 2, 2), 3);
825  FTensor::Tensor2<double, 3, 3> t_b;
826 
827  t_a(i) =
828  t_row_base * t_cauchy_stress(i, j) * t_diff_rho(j) * val * drho_coef;
829 
830  for (int cc = 0; cc != col_nb_dofs / 3; cc++) {
831 
832  t_mat(i, j) += t_a(i) * (t_invH(l, j) * t_col_diff_base(l));
833 
834  ++t_mat;
835  ++t_col_diff_base;
836  ++t_col_n;
837  }
838  ++t_row_base;
839  ++t_row_diff_n;
840  }
841 
842  for (; rr != nb_row_base; ++rr) {
843  ++t_row_base;
844  ++t_row_diff_n;
845  }
846 
847  ++t_h;
848  ++t_H;
849  ++t_invH;
850  ++rho;
851  ++t_diff_rho;
852  ++t_cauchy_stress;
853  ++det_H;
854  ++t_initial_singular_displacement;
855  }
856 
857  int *indices_ptr = &row_data.getIndices()[0];
858  if (!forcesOnlyOnEntitiesRow.empty()) {
859  iNdices.resize(row_nb_dofs);
860  noalias(iNdices) = row_data.getIndices();
861  indices_ptr = &iNdices[0];
862  auto &dofs = row_data.getFieldDofs();
863  auto dit = dofs.begin();
864  for (int ii = 0; dit != dofs.end(); dit++, ii++) {
865  if(auto dof = (*dit)) {
866  if (forcesOnlyOnEntitiesRow.find(dof->getEnt()) ==
867  forcesOnlyOnEntitiesRow.end()) {
868  iNdices[ii] = -1;
869  }
870  }
871  }
872  }
873 
874  CHKERR MatSetValues(getFEMethod()->snes_B, row_nb_dofs, indices_ptr,
875  col_nb_dofs, &*col_data.getIndices().begin(),
876  &*nA.data().begin(), ADD_VALUES);
877 
878  MoFEMFunctionReturn(0);
879 }
880 
881 OpAnalyticalSpatialTraction::OpAnalyticalSpatialTraction(
882  MoFEM::Interface &m_field, const bool &set_singular_coordinates,
883  const Range &crack_front_nodes, const Range &crack_front_nodes_edges,
884  const Range &crack_front_elements, PetscBool add_singularity,
885  const std::string field_name, const Range tris,
886  boost::ptr_vector<MethodForForceScaling> &methods_op,
887  boost::ptr_vector<NeumannForcesSurface::MethodForAnalyticalForce>
888  &analytical_force_op)
889  : FaceElementForcesAndSourcesCore::UserDataOperator(
890  field_name, UserDataOperator::OPROW),
891  tRis(tris), methodsOp(methods_op), analyticalForceOp(analytical_force_op),
892  volSideFe(m_field, set_singular_coordinates, crack_front_nodes,
893  crack_front_nodes_edges, crack_front_elements, add_singularity),
894  setSingularCoordinates(set_singular_coordinates) {}
895 
896 MoFEMErrorCode
897 OpAnalyticalSpatialTraction::doWork(int side, EntityType type,
898  DataForcesAndSourcesCore::EntData &data) {
899  MoFEMFunctionBeginHot;
900 
901  if (data.getIndices().size() == 0)
902  MoFEMFunctionReturnHot(0);
903  EntityHandle ent = getNumeredEntFiniteElementPtr()->getEnt();
904  if (tRis.find(ent) == tRis.end())
905  MoFEMFunctionReturnHot(0);
906 
907  if (type == MBVERTEX) {
908  CHKERR loopSideVolumes("ELASTIC", volSideFe);
909  }
910 
911  int rank = data.getFieldDofs()[0]->getNbOfCoeffs();
912  int nb_row_dofs = data.getIndices().size() / rank;
913 
914  Nf.resize(data.getIndices().size(), false);
915  Nf.clear();
916 
917  VectorDouble3 coords(3);
918  VectorDouble3 normal(3);
919  VectorDouble3 force(3);
920  VectorDouble3 FTforce(3);
921 
922  FTensor::Index<'i', 3> i;
923  FTensor::Index<'j', 3> j;
924  FTensor::Index<'k', 3> k;
925 
926  // If singular element apply transformation to h tensor
927  if (type == MBVERTEX && volSideFe.singularElement) {
928  double *inv_jac_ptr = &*volSideFe.invSJac.data().begin();
929  FTensor::Tensor2<double *, 3, 3> t_inv_s_jac(
930  inv_jac_ptr, &inv_jac_ptr[1], &inv_jac_ptr[2], &inv_jac_ptr[3],
931  &inv_jac_ptr[4], &inv_jac_ptr[5], &inv_jac_ptr[6], &inv_jac_ptr[7],
932  &inv_jac_ptr[8], 9);
933  FTensor::Tensor1<double *, 3> t_normal = getFTensor1Normal();
934  double nrm = 2 * getArea();
935  FTensor::Tensor1<double, 3> t_s_normal;
936  sNrm.resize(data.getN().size1(), false);
937  for (unsigned int gg = 0; gg != data.getN().size1(); gg++) {
938  // Push normal to tranformed singular element
939  double det = volSideFe.detS[gg];
940  t_s_normal(i) = det * t_inv_s_jac(j, i) * t_normal(j) / nrm;
941  sNrm[gg] = sqrt(t_s_normal(i) * t_s_normal(i));
942  ++t_inv_s_jac;
943  }
944  }
945 
946  FTensor::Tensor1<double *, 3> t_force(&force[0], &force[1], &force[2]);
947 
948  const auto &coords_at_gauss_pts = getCoordsAtGaussPts();
949  const auto &normals_at_gauss_pts = getNormalsAtGaussPts();
950 
951  for (unsigned int gg = 0; gg < data.getN().size1(); gg++) {
952 
953  // get integration weight and Jacobian of integration point (area of face)
954  double val = getGaussPts()(2, gg);
955  val *= 0.5 * cblas_dnrm2(3, &getNormalsAtGaussPts()(gg, 0), 1);
956  for (int dd = 0; dd != 3; dd++) {
957  coords[dd] = coords_at_gauss_pts(gg, dd);
958  normal[dd] = normals_at_gauss_pts(gg, dd);
959  }
960  if (volSideFe.singularElement) {
961  for (int dd = 0; dd != 3; dd++) {
962  coords[dd] += volSideFe.singularDisp(gg, dd);
963  }
964  val *= sNrm[gg];
965  }
966 
967  for (boost::ptr_vector<
968  NeumannForcesSurface::MethodForAnalyticalForce>::iterator vit =
969  analyticalForceOp.begin();
970  vit != analyticalForceOp.end(); vit++) {
971 
972  // Calculate force
973  CHKERR vit->getForce(ent, coords, normal, force);
974 
975  for (int rr = 0; rr != 3; rr++) {
976  cblas_daxpy(nb_row_dofs, val * force[rr], &data.getN()(gg, 0), 1,
977  &Nf[rr], rank);
978  }
979  }
980  }
981 
982  // Scale force using user defined scaling operator
983  CHKERR MethodForForceScaling::applyScale(getFEMethod(), methodsOp, Nf);
984 
985  // Assemble force into vector
986  CHKERR VecSetValues(getFEMethod()->snes_f, data.getIndices().size(),
987  &data.getIndices()[0], &Nf[0], ADD_VALUES);
988 
989  MoFEMFunctionReturnHot(0);
990 }
991 
992 OpAnalyticalMaterialTraction::OpAnalyticalMaterialTraction(
993  MoFEM::Interface &m_field, const bool &set_singular_coordinates,
994  const Range &crack_front_nodes, const Range &crack_front_nodes_edges,
995  const Range &crack_front_elements, PetscBool add_singularity,
996  const std::string field_name, const Range tris,
997  boost::ptr_vector<MethodForForceScaling> &methods_op,
998  boost::ptr_vector<NeumannForcesSurface::MethodForAnalyticalForce>
999  &analytical_force_op,
1000  Range *forces_only_on_entities_row)
1001  : FaceElementForcesAndSourcesCore::UserDataOperator(
1002  field_name, UserDataOperator::OPROW),
1003  tRis(tris), methodsOp(methods_op), analyticalForceOp(analytical_force_op),
1004  volSideFe(m_field, set_singular_coordinates, crack_front_nodes,
1005  crack_front_nodes_edges, crack_front_elements, add_singularity),
1006  setSingularCoordinates(set_singular_coordinates) {
1007  hG = boost::shared_ptr<MatrixDouble>(new MatrixDouble());
1008  HG = boost::shared_ptr<MatrixDouble>(new MatrixDouble());
1009  volSideFe.getOpPtrVector().push_back(
1010  new OpCalculateVectorFieldGradient<3, 3>("SPATIAL_POSITION", hG));
1011  volSideFe.getOpPtrVector().push_back(
1012  new OpCalculateVectorFieldGradient<3, 3>("MESH_NODE_POSITIONS", HG));
1013  volSideFe.meshPositionsFieldName = "NONE";
1014  if (forces_only_on_entities_row) {
1015  forcesOnlyOnEntitiesRow = *forces_only_on_entities_row;
1016  }
1017 }
1018 
1019 MoFEMErrorCode
1020 OpAnalyticalMaterialTraction::doWork(int side, EntityType type,
1021  DataForcesAndSourcesCore::EntData &data) {
1022  MoFEMFunctionBeginHot;
1023 
1024  if (data.getIndices().size() == 0)
1025  MoFEMFunctionReturnHot(0);
1026  EntityHandle ent = getNumeredEntFiniteElementPtr()->getEnt();
1027  if (tRis.find(ent) == tRis.end())
1028  MoFEMFunctionReturnHot(0);
1029 
1030  if (type == MBVERTEX) {
1031  CHKERR loopSideVolumes("MATERIAL", volSideFe);
1032  }
1033  if (type != MBVERTEX) {
1034  if (volSideFe.singularElement) {
1035  const EntityHandle ent = data.getFieldDofs()[0]->getEnt();
1036  int side_number, sense, offset;
1037  CHKERR volSideFe.mField.get_moab().side_number(
1038  volSideFe.numeredEntFiniteElementPtr->getEnt(), ent, side_number,
1039  sense, offset);
1040  data.getN() =
1041  volSideFe.getEntData(H1, type, side_number).getN(data.getBase());
1042  }
1043  } else if (volSideFe.singularElement) {
1044  // Set value of base functions at singular integration points
1045  // FIXME: That need to be moved to separate operator such that values
1046  // of base function are appropriately evaluated
1047  for (unsigned int gg = 0; gg != volSideFe.gaussPts.size2(); gg++) {
1048  for (int dd = 0; dd != 3; dd++) {
1049  volSideFe.gaussPts(dd, gg) += volSideFe.singularRefCoords(gg, dd);
1050  }
1051  }
1052  CHKERR ShapeMBTET(
1053  &*volSideFe.getEntData(H1, MBVERTEX, 0).getN(NOBASE).data().begin(),
1054  &volSideFe.gaussPts(0, 0), &volSideFe.gaussPts(1, 0),
1055  &volSideFe.gaussPts(2, 0), volSideFe.gaussPts.size2());
1056  CHKERR volSideFe.calHierarchicalBaseFunctionsOnElement(data.getBase());
1057  map<int, int> map_face_nodes;
1058  for (unsigned int gg = 0; gg != data.getN().size1(); gg++) {
1059  for (int nn = 0; nn != 3; nn++) {
1060  data.getN()(gg, nn) =
1061  volSideFe.getEntData(H1, MBVERTEX, 0)
1062  .getN(data.getBase())(gg, volSideFe.getFaceConnMap()[nn]);
1063  }
1064  }
1065  }
1066 
1067  int rank = data.getFieldDofs()[0]->getNbOfCoeffs();
1068  int nb_row_base = data.getIndices().size() / rank;
1069 
1070  Nf.resize(data.getIndices().size(), false);
1071  Nf.clear();
1072 
1073  VectorDouble3 coords(3);
1074  VectorDouble3 normal(3);
1075  VectorDouble3 force(3);
1076  VectorDouble3 FTforce(3);
1077 
1078  FTensor::Index<'i', 3> i;
1079  FTensor::Index<'j', 3> j;
1080  FTensor::Index<'k', 3> k;
1081 
1082  // If singular element apply transformation to h tensor
1083  if (type == MBVERTEX && volSideFe.singularElement) {
1084  // cerr << volSideFe.sJac << endl;
1085  // cerr << volSideFe.invSJac << endl;
1086  double *inv_jac_ptr = &*volSideFe.invSJac.data().begin();
1087  FTensor::Tensor2<double *, 3, 3> t_inv_s_jac(
1088  inv_jac_ptr, &inv_jac_ptr[1], &inv_jac_ptr[2], &inv_jac_ptr[3],
1089  &inv_jac_ptr[4], &inv_jac_ptr[5], &inv_jac_ptr[6], &inv_jac_ptr[7],
1090  &inv_jac_ptr[8], 9);
1091  auto t_hg = getFTensor2FromMat<3, 3>(*hG);
1092  auto t_normal = getFTensor1Normal();
1093  double nrm = 2 * getArea();
1094  FTensor::Tensor1<double, 3> t_s_normal;
1095  sNrm.resize(data.getN().size1(), false);
1096  for (unsigned int gg = 0; gg != data.getN().size1(); gg++) {
1097  // Push normal to transformed singular element
1098  double det = volSideFe.detS[gg];
1099  t_s_normal(i) = det * t_inv_s_jac(j, i) * t_normal(j) / nrm;
1100  sNrm[gg] = sqrt(t_s_normal(i) * t_s_normal(i));
1101  FTensor::Tensor2<double, 3, 3> t_tmp;
1102  t_tmp(i, j) = t_hg(i, k) * t_inv_s_jac(k, j);
1103  t_hg(i, j) = t_tmp(i, j);
1104  ++t_inv_s_jac;
1105  ++t_hg;
1106  }
1107  }
1108 
1109  FTensor::Tensor1<double *, 3> t_force(&force[0], &force[1], &force[2]);
1110  FTensor::Tensor1<double *, 3> t_ft_force(&FTforce[0], &FTforce[1],
1111  &FTforce[2]);
1112  auto t_hg = getFTensor2FromMat<3, 3>(*hG);
1113  auto t_Hg = getFTensor2FromMat<3, 3>(*HG);
1114  FTensor::Tensor2<double, 3, 3> t_inv_Hg;
1115  FTensor::Tensor2<double, 3, 3> t_F;
1116 
1117  // cerr << "\n\n";
1118  for (unsigned int gg = 0; gg < data.getN().size1(); gg++) {
1119 
1120  // get integration weight and Jacobian of integration point (area of face)
1121  double val = getGaussPts()(2, gg);
1122  val *= 0.5 * cblas_dnrm2(3, &getNormalsAtGaussPts()(gg, 0), 1);
1123  for (int dd = 0; dd != 3; dd++) {
1124  coords[dd] = getCoordsAtGaussPts()(gg, dd);
1125  normal[dd] = getNormalsAtGaussPts()(gg, dd);
1126  }
1127  if (volSideFe.singularElement) {
1128  for (int dd = 0; dd != 3; dd++) {
1129  coords[dd] += volSideFe.singularDisp(gg, dd);
1130  }
1131  val *= sNrm[gg];
1132  }
1133 
1134  for (boost::ptr_vector<
1135  NeumannForcesSurface::MethodForAnalyticalForce>::iterator vit =
1136  analyticalForceOp.begin();
1137  vit != analyticalForceOp.end(); vit++) {
1138 
1139  // Calculate force
1140  CHKERR vit->getForce(ent, coords, normal, force);
1141  // Invert matrix H
1142  double det;
1143  CHKERR determinantTensor3by3(t_Hg, det);
1144  CHKERR invertTensor3by3(t_Hg, det, t_inv_Hg);
1145 
1146  // Calculate gradient of deformation
1147  t_F(i, k) = t_hg(i, j) * t_inv_Hg(j, k);
1148  t_ft_force(i) = -t_F(j, i) * t_force(j);
1149 
1150  for (int rr = 0; rr != 3; rr++) {
1151  cblas_daxpy(nb_row_base, val * FTforce[rr], &data.getN()(gg, 0), 1,
1152  &Nf[rr], rank);
1153  }
1154  }
1155 
1156  ++t_hg;
1157  ++t_Hg;
1158  }
1159 
1160  // Scale force using user defined scaling operator
1161  CHKERR MethodForForceScaling::applyScale(getFEMethod(), methodsOp, Nf);
1162 
1163  // Assemble force into vector
1164  int nb_dofs = data.getIndices().size();
1165  int *indices_ptr = &data.getIndices()[0];
1166  if (!forcesOnlyOnEntitiesRow.empty()) {
1167  iNdices.resize(nb_dofs, false);
1168  noalias(iNdices) = data.getIndices();
1169  indices_ptr = &iNdices[0];
1170  auto &dofs = data.getFieldDofs();
1171  auto dit = dofs.begin();
1172  for (int ii = 0; dit != dofs.end(); dit++, ii++) {
1173  if (auto dof = (*dit)) {
1174  if (forcesOnlyOnEntitiesRow.find(dof->getEnt()) ==
1175  forcesOnlyOnEntitiesRow.end()) {
1176  iNdices[ii] = -1;
1177  }
1178  }
1179  }
1180  }
1181  CHKERR VecSetValues(getFEMethod()->snes_f, nb_dofs, indices_ptr, &Nf[0],
1182  ADD_VALUES);
1183  MoFEMFunctionReturnHot(0);
1184 }
1185 
1186 OpAleLhsWithDensitySingularElement_dX_dX::
1187  OpAleLhsWithDensitySingularElement_dX_dX(
1188  const std::string row_field, const std::string col_field,
1189  boost::shared_ptr<HookeElement::DataAtIntegrationPts> &data_at_pts,
1190  boost::shared_ptr<VectorDouble> rho_at_gauss_pts,
1191  boost::shared_ptr<MatrixDouble> rho_grad_at_gauss_pts,
1192  const double rho_n, const double rho_0,
1193  boost::shared_ptr<MatrixDouble> singular_displacement)
1194  : UserDataOperator(row_field, col_field, OPROWCOL, true),
1195  dataAtPts(data_at_pts), rhoAtGaussPtsPtr(rho_at_gauss_pts),
1196  rhoGradAtGaussPtsPtr(rho_grad_at_gauss_pts),
1197  singularDisplacement(singular_displacement), rhoN(rho_n), rHo0(rho_0) {}
1198 
1199 MoFEMErrorCode OpAleLhsWithDensitySingularElement_dX_dX::doWork(
1200  int row_side, int col_side, EntityType row_type, EntityType col_type,
1201  DataForcesAndSourcesCore::EntData &row_data,
1202  DataForcesAndSourcesCore::EntData &col_data) {
1203 
1204  MoFEMFunctionBegin;
1205  nbRows = row_data.getIndices().size();
1206  if (!nbRows)
1207  MoFEMFunctionReturnHot(0);
1208 
1209  nbCols = col_data.getIndices().size();
1210  if (!nbCols)
1211  MoFEMFunctionReturnHot(0);
1212 
1213  K.resize(nbRows, nbCols, false);
1214  K.clear();
1215 
1216  nbIntegrationPts = getGaussPts().size2();
1217  if (row_side == col_side && row_type == col_type) {
1218  isDiag = true;
1219  } else {
1220  isDiag = false;
1221  }
1222 
1223  CHKERR iNtegrate(row_data, col_data);
1224  CHKERR aSsemble(row_data, col_data);
1225 
1226  MoFEMFunctionReturn(0);
1227 }
1228 
1229 MoFEMErrorCode OpAleLhsWithDensitySingularElement_dX_dX::iNtegrate(
1230  DataForcesAndSourcesCore::EntData &row_data,
1231  DataForcesAndSourcesCore::EntData &col_data) {
1232  MoFEMFunctionBegin;
1233  auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
1234  return FTensor::Tensor2<FTensor::PackPtr<double *, 3>, 3, 3>(
1235  &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2), &m(r + 1, c + 0),
1236  &m(r + 1, c + 1), &m(r + 1, c + 2), &m(r + 2, c + 0), &m(r + 2, c + 1),
1237  &m(r + 2, c + 2));
1238  };
1239 
1240  FTensor::Index<'i', 3> i;
1241  FTensor::Index<'j', 3> j;
1242  FTensor::Index<'k', 3> k;
1243  FTensor::Index<'l', 3> l;
1244 
1245  double vol = getVolume();
1246 
1247  auto t_w = getFTensor0IntegrationWeight();
1248 
1249  auto t_row_diff_base = row_data.getFTensor1DiffN<3>();
1250  const int row_nb_base_fun = row_data.getN().size2();
1251 
1252  auto rho = getFTensor0FromVec(*rhoAtGaussPtsPtr);
1253  auto t_grad_rho = getFTensor1FromMat<3>(*rhoGradAtGaussPtsPtr);
1254 
1255  auto t_eshelby_stress =
1256  getFTensor2FromMat<3, 3>(*dataAtPts->eshelbyStressMat);
1257  auto t_invH = getFTensor2FromMat<3, 3>(*dataAtPts->invHMat);
1258  auto &det_H = *dataAtPts->detHVec;
1259 
1260  auto t_initial_singular_displacement =
1261  getFTensor1FromMat<3>(*singularDisplacement);
1262 
1263  for (int gg = 0; gg != nbIntegrationPts; ++gg) {
1264 
1265  double a = t_w * vol * det_H[gg];
1266  const double stress_dho_coef = (rhoN / rho);
1267 
1268  int rr = 0;
1269  for (; rr != nbRows / 3; ++rr) {
1270 
1271  auto t_m = get_tensor2(K, 3 * rr, 0);
1272 
1273  FTensor::Tensor1<double, 3> t_row_diff_base_pulled;
1274  t_row_diff_base_pulled(i) = t_row_diff_base(j) * t_invH(j, i);
1275 
1276  FTensor::Tensor1<double, 3> t_row_stress;
1277  t_row_stress(i) = a * t_row_diff_base_pulled(j) * t_eshelby_stress(i, j);
1278  FTensor::Tensor2<double, 3, 3> t_a;
1279  t_a(i, k) = t_row_stress(i) * stress_dho_coef * t_grad_rho(k);
1280 
1281  auto t_col_diff_base = col_data.getFTensor1DiffN<3>(gg, 0);
1282  auto t_col_base = col_data.getFTensor0N(gg, 0);
1283  for (int cc = 0; cc != nbCols / 3; ++cc) {
1284 
1285  t_m(i, k) +=
1286  t_a(i, k) * t_col_diff_base(l) * t_initial_singular_displacement(l);
1287 
1288  ++t_col_base;
1289  ++t_col_diff_base;
1290  ++t_m;
1291  }
1292  ++t_row_diff_base;
1293  }
1294 
1295  for (; rr != row_nb_base_fun; ++rr)
1296  ++t_row_diff_base;
1297 
1298  ++t_initial_singular_displacement;
1299  ++t_w;
1300  ++t_eshelby_stress;
1301  ++t_invH;
1302  ++rho;
1303  ++t_grad_rho;
1304  }
1305 
1306  MoFEMFunctionReturn(0);
1307 }
1308 
1309 MoFEMErrorCode OpAleLhsWithDensitySingularElement_dX_dX::aSsemble(
1310  DataForcesAndSourcesCore::EntData &row_data,
1311  DataForcesAndSourcesCore::EntData &col_data) {
1312  MoFEMFunctionBegin;
1313 
1314  const int *row_indices = &*row_data.getIndices().data().begin();
1315  const int *col_indices = &*col_data.getIndices().data().begin();
1316 
1317  auto &data = *dataAtPts;
1318  Mat B = getFEMethod()->ksp_B != PETSC_NULL ? getFEMethod()->ksp_B
1319  : getFEMethod()->snes_B;
1320  // assemble local matrix
1321  CHKERR MatSetValues(B, nbRows, row_indices, nbCols, col_indices,
1322  &*K.data().begin(), ADD_VALUES);
1323 
1324  if (!isDiag && sYmm) {
1325  // if not diagonal term and since global matrix is symmetric assemble
1326  // transpose term.
1327  transK.resize(K.size2(), K.size1(), false);
1328  noalias(transK) = trans(K);
1329  CHKERR MatSetValues(B, nbCols, col_indices, nbRows, row_indices,
1330  &*transK.data().begin(), ADD_VALUES);
1331  }
1332  MoFEMFunctionReturn(0);
1333 }
1334 
1335 OpAleLhsWithDensitySingularElement_dx_dX::
1336  OpAleLhsWithDensitySingularElement_dx_dX(
1337  const std::string row_field, const std::string col_field,
1338  boost::shared_ptr<HookeElement::DataAtIntegrationPts> &data_at_pts,
1339  boost::shared_ptr<VectorDouble> rho_at_gauss_pts,
1340  boost::shared_ptr<MatrixDouble> rho_grad_at_gauss_pts,
1341  const double rho_n, const double rho_0,
1342  boost::shared_ptr<MatrixDouble> singular_displacement)
1343  : UserDataOperator(row_field, col_field, OPROWCOL, false),
1344  dataAtPts(data_at_pts), rhoAtGaussPtsPtr(rho_at_gauss_pts),
1345  rhoGradAtGaussPtsPtr(rho_grad_at_gauss_pts),
1346  singularDisplacement(singular_displacement), rhoN(rho_n), rHo0(rho_0) {}
1347 
1348 MoFEMErrorCode OpAleLhsWithDensitySingularElement_dx_dX::doWork(
1349  int row_side, int col_side, EntityType row_type, EntityType col_type,
1350  DataForcesAndSourcesCore::EntData &row_data,
1351  DataForcesAndSourcesCore::EntData &col_data) {
1352 
1353  MoFEMFunctionBegin;
1354  nbRows = row_data.getIndices().size();
1355  if (!nbRows)
1356  MoFEMFunctionReturnHot(0);
1357 
1358  nbCols = col_data.getIndices().size();
1359  if (!nbCols)
1360  MoFEMFunctionReturnHot(0);
1361 
1362  K.resize(nbRows, nbCols, false);
1363  K.clear();
1364 
1365  nbIntegrationPts = getGaussPts().size2();
1366  if (row_side == col_side && row_type == col_type) {
1367  isDiag = true;
1368  } else {
1369  isDiag = false;
1370  }
1371 
1372  CHKERR iNtegrate(row_data, col_data);
1373  CHKERR aSsemble(row_data, col_data);
1374 
1375  MoFEMFunctionReturn(0);
1376 }
1377 
1378 MoFEMErrorCode OpAleLhsWithDensitySingularElement_dx_dX::iNtegrate(
1379  DataForcesAndSourcesCore::EntData &row_data,
1380  DataForcesAndSourcesCore::EntData &col_data) {
1381  MoFEMFunctionBegin;
1382  auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
1383  return FTensor::Tensor2<FTensor::PackPtr<double *, 3>, 3, 3>(
1384  &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2), &m(r + 1, c + 0),
1385  &m(r + 1, c + 1), &m(r + 1, c + 2), &m(r + 2, c + 0), &m(r + 2, c + 1),
1386  &m(r + 2, c + 2));
1387  };
1388 
1389  FTensor::Index<'i', 3> i;
1390  FTensor::Index<'j', 3> j;
1391  FTensor::Index<'k', 3> k;
1392  FTensor::Index<'l', 3> l;
1393 
1394  double vol = getVolume();
1395 
1396  auto t_w = getFTensor0IntegrationWeight();
1397 
1398  auto t_row_diff_base = row_data.getFTensor1DiffN<3>();
1399  const int row_nb_base_fun = row_data.getN().size2();
1400 
1401  auto rho = getFTensor0FromVec(*rhoAtGaussPtsPtr);
1402  auto t_grad_rho = getFTensor1FromMat<3>(*rhoGradAtGaussPtsPtr);
1403 
1404  auto t_cauchy_stress =
1405  getFTensor2SymmetricFromMat<3>(*(dataAtPts->cauchyStressMat));
1406  auto t_invH = getFTensor2FromMat<3, 3>(*dataAtPts->invHMat);
1407  auto &det_H = *dataAtPts->detHVec;
1408 
1409  auto t_initial_singular_displacement =
1410  getFTensor1FromMat<3>(*singularDisplacement);
1411 
1412  for (int gg = 0; gg != nbIntegrationPts; ++gg) {
1413 
1414  double a = t_w * vol * det_H[gg];
1415  const double stress_dho_coef = (rhoN / rho);
1416 
1417  int rr = 0;
1418  for (; rr != nbRows / 3; ++rr) {
1419 
1420  auto t_m = get_tensor2(K, 3 * rr, 0);
1421 
1422  FTensor::Tensor1<double, 3> t_row_diff_base_pulled;
1423  t_row_diff_base_pulled(i) = t_row_diff_base(j) * t_invH(j, i);
1424 
1425  FTensor::Tensor1<double, 3> t_row_stress;
1426  t_row_stress(i) = a * t_row_diff_base_pulled(j) * t_cauchy_stress(i, j);
1427  FTensor::Tensor2<double, 3, 3> t_a;
1428  t_a(i, k) = t_row_stress(i) * stress_dho_coef * t_grad_rho(k);
1429 
1430  auto t_col_diff_base = col_data.getFTensor1DiffN<3>(gg, 0);
1431  auto t_col_base = col_data.getFTensor0N(gg, 0);
1432  for (int cc = 0; cc != nbCols / 3; ++cc) {
1433 
1434  t_m(i, k) +=
1435  t_a(i, k) * t_col_diff_base(l) * t_initial_singular_displacement(l);
1436 
1437  ++t_col_base;
1438  ++t_col_diff_base;
1439  ++t_m;
1440  }
1441  ++t_row_diff_base;
1442  }
1443 
1444  for (; rr != row_nb_base_fun; ++rr)
1445  ++t_row_diff_base;
1446 
1447  ++t_initial_singular_displacement;
1448  ++t_w;
1449  ++t_cauchy_stress;
1450  ++t_invH;
1451  ++rho;
1452  ++t_grad_rho;
1453  }
1454 
1455  MoFEMFunctionReturn(0);
1456 }
1457 
1458 MoFEMErrorCode OpAleLhsWithDensitySingularElement_dx_dX::aSsemble(
1459  DataForcesAndSourcesCore::EntData &row_data,
1460  DataForcesAndSourcesCore::EntData &col_data) {
1461  MoFEMFunctionBegin;
1462 
1463  const int *row_indices = &*row_data.getIndices().data().begin();
1464  const int *col_indices = &*col_data.getIndices().data().begin();
1465 
1466  auto &data = *dataAtPts;
1467  Mat B = getFEMethod()->ksp_B != PETSC_NULL ? getFEMethod()->ksp_B
1468  : getFEMethod()->snes_B;
1469  // assemble local matrix
1470  CHKERR MatSetValues(B, nbRows, row_indices, nbCols, col_indices,
1471  &*K.data().begin(), ADD_VALUES);
1472 
1473  if (!isDiag && sYmm) {
1474  // if not diagonal term and since global matrix is symmetric assemble
1475  // transpose term.
1476  transK.resize(K.size2(), K.size1(), false);
1477  noalias(transK) = trans(K);
1478  CHKERR MatSetValues(B, nbCols, col_indices, nbRows, row_indices,
1479  &*transK.data().begin(), ADD_VALUES);
1480  }
1481  MoFEMFunctionReturn(0);
1482 }
1483 
1484 MoFEMErrorCode
1485 OpGetDensityFieldForTesting::doWork(int row_side, EntityType row_type,
1486  HookeElement::EntData &row_data) {
1487  MoFEMFunctionBegin;
1488  if (row_type != MBVERTEX)
1489  MoFEMFunctionReturnHot(0);
1490  // get number of integration points
1491  const int nb_integration_pts = getGaussPts().size2();
1492  rhoAtGaussPtsPtr->resize(nb_integration_pts, false);
1493  rhoAtGaussPtsPtr->clear();
1494  rhoGradAtGaussPtsPtr->resize(3, nb_integration_pts, false);
1495  rhoGradAtGaussPtsPtr->clear();
1496  rhoGradGradAtGaussPtsPtr->resize(9, nb_integration_pts, false);
1497  rhoGradGradAtGaussPtsPtr->clear();
1498 
1499  auto t_rho = getFTensor0FromVec(*rhoAtGaussPtsPtr);
1500  auto t_grad_rho = getFTensor1FromMat<3>(*rhoGradAtGaussPtsPtr);
1501  auto t_grad_grad_rho = getFTensor2FromMat<3, 3>(*rhoGradGradAtGaussPtsPtr);
1502 
1503  MatrixDouble &coords = *matCoordsPtr;
1504 
1505  FTensor::Index<'i', 3> i;
1506  FTensor::Index<'j', 3> j;
1507 
1508  singularDisp.resize(3, nb_integration_pts, false);
1509  singularDisp.clear();
1510 
1511  auto t_material_positions = getFTensor1FromMat<3>(coords);
1512  auto t_mwls_material_positions = getFTensor1FromMat<3>(coords);
1513  if (feSingularPtr->singularElement) {
1514  MatrixDouble singular_current_displacement;
1515 
1516  singular_current_displacement.resize(3, nb_integration_pts, false);
1517  singular_current_displacement.clear();
1518 
1519  FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_singular_displacement(
1520  &feSingularPtr->singularDisp(0, 0), &feSingularPtr->singularDisp(0, 1),
1521  &feSingularPtr->singularDisp(0, 2));
1522  auto t_inital_singular_displacement = getFTensor1FromMat<3>(singularDisp);
1523  auto t_current_singular_displacement =
1524  getFTensor1FromMat<3>(singular_current_displacement);
1525  if (!matGradPosAtPtsPtr)
1526  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1527  "Matrix for gradient of positions not allocated");
1528  auto t_H = getFTensor2FromMat<3, 3>(*matGradPosAtPtsPtr);
1529  for (int gg = 0; gg != nb_integration_pts; gg++) {
1530 
1531  t_inital_singular_displacement(i) = t_singular_displacement(i);
1532  t_current_singular_displacement(i) =
1533  t_H(i, j) * t_singular_displacement(j);
1534 
1535  t_material_positions(i) += t_current_singular_displacement(i);
1536 
1537  ++t_material_positions;
1538  ++t_H;
1539  ++t_singular_displacement;
1540  ++t_inital_singular_displacement;
1541  ++t_current_singular_displacement;
1542  }
1543  }
1544  for (int gg = 0; gg != nb_integration_pts; ++gg) {
1545  t_rho =
1546  1 + t_mwls_material_positions(i) *
1547  t_mwls_material_positions(i); // density is equal to
1548  // 1+x^2+y^2+z^2 (x,y,z coordines)
1549  t_grad_rho(i) = 2 * t_mwls_material_positions(i);
1550  t_grad_grad_rho(0, 0) = 2.;
1551  t_grad_grad_rho(1, 1) = 2.;
1552  t_grad_grad_rho(2, 2) = 2.;
1553 
1554  ++t_rho;
1555  ++t_mwls_material_positions;
1556  ++t_grad_rho;
1557  ++t_grad_grad_rho;
1558  }
1559 
1560  MoFEMFunctionReturn(0);
1561 }
1562 
1563 MoFEMErrorCode
1564 OpSetTagRangeOnSkin::doWork(int side, EntityType type,
1565  DataForcesAndSourcesCore::EntData &data) {
1566  MoFEMFunctionBeginHot;
1567 
1568  if (data.getFieldData().size() == 0)
1569  MoFEMFunctionReturnHot(0);
1570 
1571  int tag_val = tagValue;
1572  EntityHandle tri_ent = getNumeredEntFiniteElementPtr()->getEnt();
1573  if (rangeToTag.find(tri_ent) == rangeToTag.end())
1574  tag_val = 0;
1575 
1576  int def_VAL = 0;
1577  Tag th;
1578  CHKERR postProcMesh.tag_get_handle(tagName.c_str(), 1, MB_TYPE_INTEGER, th,
1579  MB_TAG_CREAT | MB_TAG_SPARSE, &def_VAL);
1580 
1581  const int nb_gauss_pts = mapGaussPts.size();
1582  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
1583 
1584  // check for existing mark
1585  int check_data;
1586  CHKERR postProcMesh.tag_get_data(th, &mapGaussPts[gg], 1, &check_data);
1587 
1588  if (check_data == 0)
1589  CHKERR postProcMesh.tag_set_data(th, &mapGaussPts[gg], 1, &tag_val);
1590  }
1591  MoFEMFunctionReturnHot(0);
1592 }
1593 
1594 } // namespace FractureMechanics
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
FractureMechanics::OpSetTagRangeOnSkin::doWork
MoFEMErrorCode doWork(int row_side, EntityType row_type, DataForcesAndSourcesCore::EntData &row_data)
Definition: CrackFrontElement.cpp:1564
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition: UnknownInterface.hpp:93
MoFEM::K
VectorDouble K
Definition: Projection10NodeCoordsOnField.cpp:125
FractureMechanics::OpAnalyticalMaterialTraction::forcesOnlyOnEntitiesRow
Range forcesOnlyOnEntitiesRow
Definition: CrackFrontElement.hpp:458
FractureMechanics::OpRhsBoneExplicitDerivariveWithHooke::commonData
HookeElement::DataAtIntegrationPts & commonData
Definition: CrackFrontElement.hpp:579
FractureMechanics::OpLhsBoneExplicitDerivariveWithHooke_dX::nA
MatrixDouble nA
Definition: CrackFrontElement.hpp:624
FractureMechanics::OpAnalyticalSpatialTraction::sNrm
VectorDouble sNrm
Length of the normal vector.
Definition: CrackFrontElement.hpp:417
FractureMechanics::OpGetDensityFieldForTesting::rhoAtGaussPtsPtr
boost::shared_ptr< VectorDouble > rhoAtGaussPtsPtr
Definition: CrackFrontElement.hpp:763
H1
@ H1
continuous field
Definition: definitions.h:85
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:1644
FractureMechanics::OpAnalyticalMaterialTraction::hG
boost::shared_ptr< MatrixDouble > hG
spatial deformation gradient
Definition: CrackFrontElement.hpp:454
FractureMechanics::OpLhsBoneExplicitDerivariveWithHooke_dx::singularDispl
MatrixDouble & singularDispl
Definition: CrackFrontElement.hpp:644
MoFEM::Types::VectorDouble3
VectorBoundedArray< double, 3 > VectorDouble3
Definition: Types.hpp:92
FTensor::Tensor1
Definition: Tensor1_value.hpp:8
EntityHandle
MoFEM::MeshRefinement::addVerticesInTheMiddleOfEdges
MoFEMErrorCode addVerticesInTheMiddleOfEdges(const EntityHandle meshset, const BitRefLevel &bit, const bool recursive=false, int verb=QUIET, EntityHandle start_v=0)
make vertices in the middle of edges in meshset and add them to refinement levels defined by bit
Definition: MeshRefinement.cpp:42
FractureMechanics::OpRhsBoneExplicitDerivariveWithHooke::iNdices
ublas::vector< int > iNdices
Definition: CrackFrontElement.hpp:586
FractureMechanics::OpAleLhsWithDensitySingularElement_dx_dX::nbRows
int nbRows
number of dofs on rows
Definition: CrackFrontElement.hpp:730
FractureMechanics::CrackFrontSingularBase::singularElement
bool singularElement
Definition: CrackFrontElement.hpp:47
rho
double rho
Definition: plastic.cpp:140
MoFEM::OpCalculateVectorFieldGradient_General< Tensor_Dim0, Tensor_Dim1, double, ublas::row_major, DoubleAllocator >::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
calculate values of vector field at integration points
Definition: UserDataOperators.hpp:1442
FractureMechanics::OpLhsBoneExplicitDerivariveWithHooke_dX::diffDiffRhoAtGaussPts
boost::shared_ptr< MatrixDouble > diffDiffRhoAtGaussPts
Definition: CrackFrontElement.hpp:617
FractureMechanics::OpLhsBoneExplicitDerivariveWithHooke_dx::forcesOnlyOnEntitiesRow
Range forcesOnlyOnEntitiesRow
Definition: CrackFrontElement.hpp:648
FractureMechanics::OpRhsBoneExplicitDerivariveWithHooke::rhoAtGaussPts
boost::shared_ptr< VectorDouble > rhoAtGaussPts
Definition: CrackFrontElement.hpp:580
NOBASE
@ NOBASE
Definition: definitions.h:59
Mortar.hpp
FractureMechanics::OpAleLhsWithDensitySingularElement_dX_dX::nbRows
int nbRows
number of dofs on rows
Definition: CrackFrontElement.hpp:684
FractureMechanics::OpRhsBoneExplicitDerivariveWithHooke::boneN
const double boneN
Definition: CrackFrontElement.hpp:584
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
FractureMechanics::OpAnalyticalSpatialTraction::analyticalForceOp
boost::ptr_vector< NeumannForcesSurface::MethodForAnalyticalForce > & analyticalForceOp
Definition: CrackFrontElement.hpp:410
FractureMechanics::OpLhsBoneExplicitDerivariveWithHooke_dX::rhoAtGaussPts
boost::shared_ptr< VectorDouble > rhoAtGaussPts
Definition: CrackFrontElement.hpp:615
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
FractureMechanics::OpAnalyticalMaterialTraction::analyticalForceOp
boost::ptr_vector< NeumannForcesSurface::MethodForAnalyticalForce > & analyticalForceOp
Definition: CrackFrontElement.hpp:450
FractureMechanics::OpAleLhsWithDensitySingularElement_dx_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: CrackFrontElement.cpp:1348
FractureMechanics::OpAnalyticalSpatialTraction::methodsOp
boost::ptr_vector< MethodForForceScaling > & methodsOp
Definition: CrackFrontElement.hpp:408
FractureMechanics::OpGetDensityFieldForTesting::doWork
MoFEMErrorCode doWork(int row_side, EntityType row_type, HookeElement::EntData &row_data)
Definition: CrackFrontElement.cpp:1485
FractureMechanics::OpAleLhsWithDensitySingularElement_dx_dX::iNtegrate
MoFEMErrorCode iNtegrate(DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
Definition: CrackFrontElement.cpp:1378
FractureMechanics::OpLhsBoneExplicitDerivariveWithHooke_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: CrackFrontElement.cpp:765
FractureMechanics::OpLhsBoneExplicitDerivariveWithHooke_dx::nA
MatrixDouble nA
Definition: CrackFrontElement.hpp:650
MoFEM::th
Tag th
Definition: Projection10NodeCoordsOnField.cpp:122
FractureMechanics::OpAleLhsWithDensitySingularElement_dX_dX::aSsemble
MoFEMErrorCode aSsemble(DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
Definition: CrackFrontElement.cpp:1309
QUAD_::order
int order
Definition: quad.h:28
FractureMechanics::OpAleLhsWithDensitySingularElement_dX_dX::rhoAtGaussPtsPtr
boost::shared_ptr< VectorDouble > rhoAtGaussPtsPtr
Definition: CrackFrontElement.hpp:667
MoFEM.hpp
FractureMechanics::OpAnalyticalMaterialTraction::volSideFe
CrackFrontSingularBase< FirendVolumeOnSide, VolumeElementForcesAndSourcesCore > volSideFe
Definition: CrackFrontElement.hpp:453
BasicFiniteElements.hpp
J
FTensor::Index< 'J', DIM1 > J
Definition: level_set.cpp:30
MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator::getVolume
double getVolume() const
element volume (linear geometry)
Definition: VolumeElementForcesAndSourcesCore.hpp:161
FractureMechanics::OpLhsBoneExplicitDerivariveWithHooke_dx::rhoAtGaussPts
boost::shared_ptr< VectorDouble > rhoAtGaussPts
Definition: CrackFrontElement.hpp:641
FractureMechanics::OpAleLhsWithDensitySingularElement_dx_dX::aSsemble
MoFEMErrorCode aSsemble(DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
Definition: CrackFrontElement.cpp:1458
FractureMechanics::OpAleLhsWithDensitySingularElement_dx_dX::dataAtPts
boost::shared_ptr< HookeElement::DataAtIntegrationPts > dataAtPts
Definition: CrackFrontElement.hpp:725
FractureMechanics::OpAnalyticalMaterialTraction::Nf
VectorDouble Nf
Definition: CrackFrontElement.hpp:472
FractureMechanics::OpAleLhsWithDensitySingularElement_dX_dX::singularDisplacement
boost::shared_ptr< MatrixDouble > singularDisplacement
Definition: CrackFrontElement.hpp:669
FractureMechanics::OpAleLhsWithDensitySingularElement_dX_dX::iNtegrate
MoFEMErrorCode iNtegrate(DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
Definition: CrackFrontElement.cpp:1229
FractureMechanics::OpLhsBoneExplicitDerivariveWithHooke_dx::diffRhoAtGaussPts
boost::shared_ptr< MatrixDouble > diffRhoAtGaussPts
Definition: CrackFrontElement.hpp:642
MoFEM::VecSetValues
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
Definition: EntitiesFieldData.hpp:1589
sdf.r
int r
Definition: sdf.py:8
FTensor::Tensor2< double *, 3, 3 >
FractureMechanics::CrackFrontSingularBase::invSJac
MatrixDouble invSJac
Definition: CrackFrontElement.hpp:42
FractureMechanics::OpSetTagRangeOnSkin::mapGaussPts
std::vector< EntityHandle > & mapGaussPts
Definition: CrackFrontElement.hpp:797
order
constexpr int order
Definition: dg_projection.cpp:18
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
FractureMechanics::OpRhsBoneExplicitDerivariveWithHooke::doWork
MoFEMErrorCode doWork(int row_side, EntityType row_type, DataForcesAndSourcesCore::EntData &row_data)
Definition: CrackFrontElement.cpp:518
FractureMechanics::OpLhsBoneExplicitDerivariveWithHooke_dx::commonData
HookeElement::DataAtIntegrationPts & commonData
Definition: CrackFrontElement.hpp:640
FractureMechanics::OpAnalyticalMaterialTraction::HG
boost::shared_ptr< MatrixDouble > HG
spatial deformation gradient
Definition: CrackFrontElement.hpp:455
I
constexpr IntegrationType I
Definition: operators_tests.cpp:31
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:1588
MoFEM::ForcesAndSourcesCore::UserDataOperator::getGaussPts
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
Definition: ForcesAndSourcesCore.hpp:1236
FractureMechanics::OpRhsBoneExplicitDerivariveWithHooke::tetsInBlock
Range tetsInBlock
Definition: CrackFrontElement.hpp:582
FractureMechanics::OpAleLhsWithDensitySingularElement_dx_dX::K
MatrixDouble K
Definition: CrackFrontElement.hpp:721
FTensor::Tensor2::data
T data[Tensor_Dim0][Tensor_Dim1]
Definition: Tensor2_value.hpp:18
MoFEM::MeshRefinement::refineTets
MoFEMErrorCode refineTets(const EntityHandle meshset, const BitRefLevel &bit, int verb=QUIET, const bool debug=false)
refine TET in the meshset
Definition: MeshRefinement.cpp:197
FractureMechanics::OpLhsBoneExplicitDerivariveWithHooke_dX::iNdices
ublas::vector< int > iNdices
Definition: CrackFrontElement.hpp:623
c
const double c
speed of light (cm/ns)
Definition: initial_diffusion.cpp:39
FractureMechanics::OpAnalyticalSpatialTraction::doWork
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Definition: CrackFrontElement.cpp:897
FractureMechanics::OpAnalyticalMaterialTraction::iNdices
VectorInt iNdices
Definition: CrackFrontElement.hpp:460
NeumannForcesSurface::MethodForAnalyticalForce
Analytical force method.
Definition: SurfacePressure.hpp:21
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
FractureMechanics::OpLhsBoneExplicitDerivariveWithHooke_dX::forcesOnlyOnEntitiesRow
Range forcesOnlyOnEntitiesRow
Definition: CrackFrontElement.hpp:622
FractureMechanics::OpAleLhsWithDensitySingularElement_dx_dX::isDiag
bool isDiag
true if this block is on diagonal
Definition: CrackFrontElement.hpp:733
FractureMechanics::OpAnalyticalMaterialTraction::sNrm
VectorDouble sNrm
Length of the normal vector.
Definition: CrackFrontElement.hpp:456
FractureMechanics::OpAnalyticalMaterialTraction::tRis
const Range tRis
Definition: CrackFrontElement.hpp:447
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
a
constexpr double a
Definition: approx_sphere.cpp:30
FractureMechanics::OpAleLhsWithDensitySingularElement_dX_dX::OpAleLhsWithDensitySingularElement_dX_dX
OpAleLhsWithDensitySingularElement_dX_dX(const std::string row_field, const std::string col_field, boost::shared_ptr< HookeElement::DataAtIntegrationPts > &data_at_pts, boost::shared_ptr< VectorDouble > rho_at_gauss_pts, boost::shared_ptr< MatrixDouble > rho_grad_at_gauss_pts, const double rho_n, const double rho_0, boost::shared_ptr< MatrixDouble > singular_displacement)
Definition: CrackFrontElement.cpp:1187
MoFEM::ForcesAndSourcesCore::UserDataOperator
Definition: ForcesAndSourcesCore.hpp:549
MoFEM::MeshRefinement
Mesh refinement interface.
Definition: MeshRefinement.hpp:26
MoFEM::CoreTmp
Definition: Core.hpp:36
Hooke.hpp
Implementation of linear elastic material.
convert.type
type
Definition: convert.py:64
FractureMechanics::OpLhsBoneExplicitDerivariveWithHooke_dx::tetsInBlock
Range tetsInBlock
Definition: CrackFrontElement.hpp:645
FractureMechanics::OpRhsBoneExplicitDerivariveWithHooke::forcesOnlyOnEntitiesRow
Range forcesOnlyOnEntitiesRow
Definition: CrackFrontElement.hpp:585
MoFEM::getFTensor0FromVec
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:135
FractureMechanics::OpAleLhsWithDensitySingularElement_dX_dX::nbCols
int nbCols
number if dof on column
Definition: CrackFrontElement.hpp:685
FractureMechanics::CrackFrontSingularBase::singularDisp
MatrixDouble singularDisp
Definition: CrackFrontElement.hpp:43
FractureMechanics::OpSetTagRangeOnSkin::postProcMesh
moab::Interface & postProcMesh
Definition: CrackFrontElement.hpp:796
FractureMechanics::OpSetTagRangeOnSkin::rangeToTag
Range rangeToTag
Definition: CrackFrontElement.hpp:799
FractureMechanics::OpRhsBoneExplicitDerivariveWithHooke::nF
VectorDouble nF
Definition: CrackFrontElement.hpp:594
FractureMechanics::OpAnalyticalSpatialTraction::volSideFe
CrackFrontSingularBase< VolumeElementForcesAndSourcesCoreOnSide, VolumeElementForcesAndSourcesCore > volSideFe
Definition: CrackFrontElement.hpp:414
FractureMechanics::CrackFrontSingularBase::detS
VectorDouble detS
Definition: CrackFrontElement.hpp:45
MoFEM::ForcesAndSourcesCore::UserDataOperator::getNumeredEntFiniteElementPtr
boost::shared_ptr< const NumeredEntFiniteElement > getNumeredEntFiniteElementPtr() const
Return raw pointer to NumeredEntFiniteElement.
Definition: ForcesAndSourcesCore.hpp:1000
QUAD_::npoints
int npoints
Definition: quad.h:29
FractureMechanics::OpSetTagRangeOnSkin::tagValue
int tagValue
Definition: CrackFrontElement.hpp:800
ShapeVolumeMBTET
double ShapeVolumeMBTET(double *diffN, const double *coords)
calculate TET volume
Definition: fem_tools.c:298
MoFEM::FaceElementForcesAndSourcesCore
Face finite element.
Definition: FaceElementForcesAndSourcesCore.hpp:23
NonlinearElasticElement::OpGetDataAtGaussPts::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
operator calculating deformation gradient
Definition: NonLinearElasticElement.cpp:95
FTensor::Tensor4
Definition: Tensor4_value.hpp:18
CrackFrontElement.hpp
FractureMechanics::OpLhsBoneExplicitDerivariveWithHooke_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: CrackFrontElement.cpp:600
FractureMechanics::OpAleLhsWithDensitySingularElement_dx_dX::rhoAtGaussPtsPtr
boost::shared_ptr< VectorDouble > rhoAtGaussPtsPtr
Definition: CrackFrontElement.hpp:713
MoFEM::VolumeElementForcesAndSourcesCore
Volume finite element base.
Definition: VolumeElementForcesAndSourcesCore.hpp:26
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
QUAD_3D_TABLE
static QUAD *const QUAD_3D_TABLE[]
Definition: quad.h:187
FractureMechanics::OpAleLhsWithDensitySingularElement_dx_dX::rhoN
const double rhoN
Definition: CrackFrontElement.hpp:717
FractureMechanics::OpAleLhsWithDensitySingularElement_dx_dX::nbCols
int nbCols
number if dof on column
Definition: CrackFrontElement.hpp:731
FractureMechanics::mapRefCoords
static std::map< long int, MatrixDouble > mapRefCoords
Definition: CrackFrontElement.cpp:53
EntData
EntitiesFieldData::EntData EntData
Definition: child_and_parent.cpp:37
field_name
constexpr auto field_name
Definition: poisson_2d_homogeneous.cpp:13
FTensor::Index
Definition: Index.hpp:23
MoFEM::ForcesAndSourcesCore::UserDataOperator::getFEMethod
const FEMethod * getFEMethod() const
Return raw pointer to Finite Element Method object.
Definition: ForcesAndSourcesCore.hpp:1042
MoFEM::OpCalculateVectorFieldGradient
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Definition: UserDataOperators.hpp:1535
convert.n
n
Definition: convert.py:82
FractureMechanics::OpAleLhsWithDensitySingularElement_dX_dX::rhoGradAtGaussPtsPtr
boost::shared_ptr< MatrixDouble > rhoGradAtGaussPtsPtr
Definition: CrackFrontElement.hpp:668
FractureMechanics::CrackFrontSingularBase
Definition: CrackFrontElement.hpp:35
FractureMechanics::OpAleLhsWithDensitySingularElement_dX_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: CrackFrontElement.cpp:1199
FractureMechanics::OpAnalyticalMaterialTraction::doWork
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Definition: CrackFrontElement.cpp:1020
FractureMechanics::OpLhsBoneExplicitDerivariveWithHooke_dX::tetsInBlock
Range tetsInBlock
Definition: CrackFrontElement.hpp:619
cholesky.hpp
cholesky decomposition
QUAD_3D_TABLE_SIZE
#define QUAD_3D_TABLE_SIZE
Definition: quad.h:186
NonlinearElasticElement::MyVolumeFE
definition of volume element
Definition: NonLinearElasticElement.hpp:31
MoFEM::determinantTensor3by3
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
Definition: Templates.hpp:1540
FractureMechanics::OpGetDensityFieldForTesting::rhoGradGradAtGaussPtsPtr
boost::shared_ptr< MatrixDouble > rhoGradGradAtGaussPtsPtr
Definition: CrackFrontElement.hpp:765
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
FractureMechanics::OpAnalyticalMaterialTraction::OpAnalyticalMaterialTraction
OpAnalyticalMaterialTraction(MoFEM::Interface &m_field, const bool &set_singular_coordinates, const Range &crack_front_nodes, const Range &crack_front_nodes_edges, const Range &crack_front_elements, PetscBool add_singularity, const std::string field_name, const Range tris, boost::ptr_vector< MethodForForceScaling > &methods_op, boost::ptr_vector< NeumannForcesSurface::MethodForAnalyticalForce > &analytical_force_op, Range *forces_only_on_entities_row=NULL)
Definition: CrackFrontElement.cpp:992
FractureMechanics::OpAleLhsWithDensitySingularElement_dX_dX::transK
MatrixDouble transK
Definition: CrackFrontElement.hpp:676
FractureMechanics::OpAleLhsWithDensitySingularElement_dx_dX::singularDisplacement
boost::shared_ptr< MatrixDouble > singularDisplacement
Definition: CrackFrontElement.hpp:715
FractureMechanics::OpAleLhsWithDensitySingularElement_dX_dX::nbIntegrationPts
int nbIntegrationPts
number of integration points
Definition: CrackFrontElement.hpp:686
FractureMechanics::OpLhsBoneExplicitDerivariveWithHooke_dx::boneN
const double boneN
Definition: CrackFrontElement.hpp:647
FractureMechanics::OpAnalyticalSpatialTraction::OpAnalyticalSpatialTraction
OpAnalyticalSpatialTraction(MoFEM::Interface &m_field, const bool &set_singular_coordinates, const Range &crack_front_nodes, const Range &crack_front_nodes_edges, const Range &crack_front_elements, PetscBool add_singularity, const std::string field_name, const Range tris, boost::ptr_vector< MethodForForceScaling > &methods_op, boost::ptr_vector< NeumannForcesSurface::MethodForAnalyticalForce > &analytical_force_op)
Definition: CrackFrontElement.cpp:881
AnalyticalFun.hpp
FractureMechanics::OpGetDensityFieldForTesting::rhoGradAtGaussPtsPtr
boost::shared_ptr< MatrixDouble > rhoGradAtGaussPtsPtr
Definition: CrackFrontElement.hpp:764
MoFEM::ForcesAndSourcesCore::UserDataOperator::getFEEntityHandle
EntityHandle getFEEntityHandle() const
Return finite element entity handle.
Definition: ForcesAndSourcesCore.hpp:1004
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1148
FractureMechanics::OpAleLhsWithDensitySingularElement_dx_dX::transK
MatrixDouble transK
Definition: CrackFrontElement.hpp:722
FractureMechanics::OpAleLhsWithDensitySingularElement_dx_dX::OpAleLhsWithDensitySingularElement_dx_dX
OpAleLhsWithDensitySingularElement_dx_dX(const std::string row_field, const std::string col_field, boost::shared_ptr< HookeElement::DataAtIntegrationPts > &data_at_pts, boost::shared_ptr< VectorDouble > rho_at_gauss_pts, boost::shared_ptr< MatrixDouble > rho_grad_at_gauss_pts, const double rho_n, const double rho_0, boost::shared_ptr< MatrixDouble > singular_displacement)
Definition: CrackFrontElement.cpp:1336
FractureMechanics::OpAnalyticalSpatialTraction::tRis
const Range tRis
Definition: CrackFrontElement.hpp:407
FractureMechanics::OpLhsBoneExplicitDerivariveWithHooke_dX::boneN
const double boneN
Definition: CrackFrontElement.hpp:621
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
FractureMechanics::OpLhsBoneExplicitDerivariveWithHooke_dX::OpLhsBoneExplicitDerivariveWithHooke_dX
OpLhsBoneExplicitDerivariveWithHooke_dX(HookeElement::DataAtIntegrationPts &common_data, boost::shared_ptr< VectorDouble > rho_at_gauss_pts, boost::shared_ptr< MatrixDouble > diff_rho_at_gauss_pts, boost::shared_ptr< MatrixDouble > diff_diff_rho_at_gauss_pts, MatrixDouble &singular_displ, Range tets_in_block, const double rho_0, const double bone_n, Range *forces_only_on_entities_row=NULL)
Definition: CrackFrontElement.cpp:581
ShapeDiffMBTET
PetscErrorCode ShapeDiffMBTET(double *diffN)
calculate derivatives of shape functions
Definition: fem_tools.c:319
FractureMechanics::OpGetDensityFieldForTesting::matCoordsPtr
boost::shared_ptr< MatrixDouble > matCoordsPtr
Definition: CrackFrontElement.hpp:762
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
ComplexConstArea.hpp
FractureMechanics::OpAnalyticalSpatialTraction::Nf
VectorDouble Nf
Definition: CrackFrontElement.hpp:429
FractureMechanics::OpAleLhsWithDensitySingularElement_dX_dX::rhoN
const double rhoN
Definition: CrackFrontElement.hpp:671
FractureMechanics::OpSetTagRangeOnSkin::tagName
const string tagName
Definition: CrackFrontElement.hpp:798
FractureMechanics::TwoType
Definition: CrackFrontElement.hpp:25
ConstantArea.hpp
FractureMechanics::OpAleLhsWithDensitySingularElement_dX_dX::isDiag
bool isDiag
true if this block is on diagonal
Definition: CrackFrontElement.hpp:687
FractureMechanics::OpAleLhsWithDensitySingularElement_dx_dX::rhoGradAtGaussPtsPtr
boost::shared_ptr< MatrixDouble > rhoGradAtGaussPtsPtr
Definition: CrackFrontElement.hpp:714
sdf_wavy_2d.w
int w
Definition: sdf_wavy_2d.py:6
m
FTensor::Index< 'm', 3 > m
Definition: shallow_wave.cpp:80
MoFEM::BitRefManager
Managing BitRefLevels.
Definition: BitRefManager.hpp:21
FractureMechanics::OpLhsBoneExplicitDerivariveWithHooke_dx::iNdices
ublas::vector< int > iNdices
Definition: CrackFrontElement.hpp:649
MoFEM::Types::BitRefLevel
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:40
FractureMechanics::OpAleLhsWithDensitySingularElement_dX_dX::dataAtPts
boost::shared_ptr< HookeElement::DataAtIntegrationPts > dataAtPts
Definition: CrackFrontElement.hpp:679
FractureMechanics::OpLhsBoneExplicitDerivariveWithHooke_dx::OpLhsBoneExplicitDerivariveWithHooke_dx
OpLhsBoneExplicitDerivariveWithHooke_dx(HookeElement::DataAtIntegrationPts &common_data, boost::shared_ptr< VectorDouble > rho_at_gauss_pts, boost::shared_ptr< MatrixDouble > diff_rho_at_gauss_pts, boost::shared_ptr< MatrixDouble > diff_diff_rho_at_gauss_pts, MatrixDouble &singular_displ, Range tets_in_block, const double rho_0, const double bone_n, Range *forces_only_on_entities_row=NULL)
Definition: CrackFrontElement.cpp:746
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:453
MethodForForceScaling::applyScale
static MoFEMErrorCode applyScale(const FEMethod *fe, boost::ptr_vector< MethodForForceScaling > &methods_op, VectorDouble &nf)
Definition: MethodForForceScaling.hpp:21
FractureMechanics::OpGetDensityFieldForTesting::singularDisp
MatrixDouble & singularDisp
Definition: CrackFrontElement.hpp:767
NeoHookean.hpp
Implementation of Neo-Hookean elastic material.
FractureMechanics::OpAleLhsWithDensitySingularElement_dX_dX::K
MatrixDouble K
Definition: CrackFrontElement.hpp:675
FractureMechanics::OpAleLhsWithDensitySingularElement_dx_dX::nbIntegrationPts
int nbIntegrationPts
number of integration points
Definition: CrackFrontElement.hpp:732
k
FTensor::Index< 'k', 3 > k
Definition: matrix_function.cpp:20
FractureMechanics::OpGetDensityFieldForTesting::matGradPosAtPtsPtr
boost::shared_ptr< MatrixDouble > matGradPosAtPtsPtr
Definition: CrackFrontElement.hpp:768
FractureMechanics::OpLhsBoneExplicitDerivariveWithHooke_dX::diffRhoAtGaussPts
boost::shared_ptr< MatrixDouble > diffRhoAtGaussPts
Definition: CrackFrontElement.hpp:616
FractureMechanics::OpLhsBoneExplicitDerivariveWithHooke_dX::singularDispl
MatrixDouble & singularDispl
Definition: CrackFrontElement.hpp:618
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
MOFEM_NOT_IMPLEMENTED
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:32
FractureMechanics::OpLhsBoneExplicitDerivariveWithHooke_dX::commonData
HookeElement::DataAtIntegrationPts & commonData
Definition: CrackFrontElement.hpp:614
FractureMechanics::OpGetDensityFieldForTesting::feSingularPtr
boost::shared_ptr< CrackFrontElement > feSingularPtr
Definition: CrackFrontElement.hpp:766
H
double H
Hardening.
Definition: plastic.cpp:124
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
l
FTensor::Index< 'l', 3 > l
Definition: matrix_function.cpp:21
FractureMechanics
Definition: AnalyticalFun.hpp:15
FractureMechanics::OpAnalyticalMaterialTraction::methodsOp
boost::ptr_vector< MethodForForceScaling > & methodsOp
Definition: CrackFrontElement.hpp:448
FractureMechanics::OpRhsBoneExplicitDerivariveWithHooke::rhoGradAtGaussPts
boost::shared_ptr< MatrixDouble > rhoGradAtGaussPts
Definition: CrackFrontElement.hpp:581
ShapeMBTET
PetscErrorCode ShapeMBTET(double *N, const double *G_X, const double *G_Y, const double *G_Z, int DIM)
calculate shape functions
Definition: fem_tools.c:306
quad.h
Generated by Doxygen 1.8.17 and hosted at University of Glasgow