v0.13.1
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>
22using 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
40extern "C" {
41#include <phg-quadrule/quad.h>
42}
43
44#include <NeoHookean.hpp>
45#include <Hooke.hpp>
46
47#include <CrackFrontElement.hpp>
48
49namespace po = boost::program_options;
50
51namespace FractureMechanics {
52
53static std::map<long int, MatrixDouble> mapRefCoords;
54
55template <>
56template <>
61 int order) {
62 return -1;
63}
64
65template <>
66template <>
71 int order) {
72
74
75 auto set_gauss_pts = [this]() {
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);
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
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)
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
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
261}
262
266
267 if (!singularElement) {
269 }
270
271 const int nb_dofs = data.getFieldData().size();
272 if (!nb_dofs && type == this->zeroType) {
273 this->dataPtr->resize(9, 0, false);
275 }
276 if (!nb_dofs) {
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);
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
320}
321
326
327 if (!singularElement) {
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) {
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>();
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>();
371 &valuesAtGaussPts[gg][1],
372 &valuesAtGaussPts[gg][2]);
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
402}
403
408 if (!singularElement || type != MBVERTEX) {
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);
419 auto t_H = getFTensor2FromMat<3, 3>(*H);
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 }
431}
432
438
439 if (!singularElement)
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();
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)
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)
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
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
502}
503
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)
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
519 int row_side, EntityType row_type,
523 const int nb_dofs = row_data.getIndices().size();
524 if (!nb_dofs)
526 if (tetsInBlock.find(getFEEntityHandle()) == tetsInBlock.end()) {
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);
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);
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()) ==
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);
578}
579
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)
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
601 int row_side, int col_side, EntityType row_type, EntityType col_type,
605 if (tetsInBlock.find(getFEEntityHandle()) == tetsInBlock.end()) {
607 }
608
609 const int row_nb_dofs = row_data.getIndices().size();
610 if (!row_nb_dofs)
612
613 const int col_nb_dofs = col_data.getIndices().size();
614 if (!col_nb_dofs)
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);
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
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);
654 t_F_dX(i, j, k, l) = -(t_h(i, m) * t_invH(m, k)) * t_H(l, j);
655
657 t_energy_dX(k, l) = t_F_dX(i, j, k, l) * t_cauchy_stress(i, j);
658
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++) {
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
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()) ==
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);
743}
744
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)
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
766 int row_side, int col_side, EntityType row_type, EntityType col_type,
770 if (tetsInBlock.find(getFEEntityHandle()) == tetsInBlock.end()) {
772 }
773 const int row_nb_dofs = row_data.getIndices().size();
774 if (!row_nb_dofs) {
776 }
777 const int col_nb_dofs = col_data.getIndices().size();
778 if (!col_nb_dofs) {
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));
799
800 auto t_diff_rho = getFTensor1FromMat<3>(*diffRhoAtGaussPts);
801
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
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);
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()) ==
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
879}
880
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)
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
900
901 if (data.getIndices().size() == 0)
903 EntityHandle ent = getNumeredEntFiniteElementPtr()->getEnt();
904 if (tRis.find(ent) == tRis.end())
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);
919 VectorDouble3 force(3);
920 VectorDouble3 FTforce(3);
921
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();
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();
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 }
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<
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
984
985 // Assemble force into vector
986 CHKERR VecSetValues(getFEMethod()->snes_f, data.getIndices().size(),
987 &data.getIndices()[0], &Nf[0], ADD_VALUES);
988
990}
991
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)
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
1023
1024 if (data.getIndices().size() == 0)
1026 EntityHandle ent = getNumeredEntFiniteElementPtr()->getEnt();
1027 if (tRis.find(ent) == tRis.end())
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 }
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);
1075 VectorDouble3 force(3);
1076 VectorDouble3 FTforce(3);
1077
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();
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));
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);
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<
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
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()) ==
1176 iNdices[ii] = -1;
1177 }
1178 }
1179 }
1180 }
1181 CHKERR VecSetValues(getFEMethod()->snes_f, nb_dofs, indices_ptr, &Nf[0],
1182 ADD_VALUES);
1184}
1185
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
1200 int row_side, int col_side, EntityType row_type, EntityType col_type,
1203
1205 nbRows = row_data.getIndices().size();
1206 if (!nbRows)
1208
1209 nbCols = col_data.getIndices().size();
1210 if (!nbCols)
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
1227}
1228
1233 auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
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
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
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);
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
1307}
1308
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 }
1333}
1334
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
1349 int row_side, int col_side, EntityType row_type, EntityType col_type,
1352
1354 nbRows = row_data.getIndices().size();
1355 if (!nbRows)
1357
1358 nbCols = col_data.getIndices().size();
1359 if (!nbCols)
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
1376}
1377
1382 auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
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
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
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);
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
1456}
1457
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 }
1482}
1483
1486 HookeElement::EntData &row_data) {
1488 if (row_type != MBVERTEX)
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
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
1561}
1562
1567
1568 if (data.getFieldData().size() == 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 }
1592}
1593
1594} // namespace FractureMechanics
static Index< 'J', 3 > J
static Index< 'K', 3 > K
Implementation of linear elastic material.
ForcesAndSourcesCore::UserDataOperator UserDataOperator
Implementation of Neo-Hookean elastic material.
constexpr double a
EntitiesFieldData::EntData EntData
cholesky decomposition
T data[Tensor_Dim0][Tensor_Dim1]
@ NOBASE
Definition: definitions.h:59
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
@ H1
continuous field
Definition: definitions.h:85
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:32
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
constexpr IntegrationType I
Definition: elastic.cpp:31
const int dim
PetscErrorCode ShapeDiffMBTET(double *diffN)
calculate derivatives of shape functions
Definition: fem_tools.c:319
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
double ShapeVolumeMBTET(double *diffN, const double *coords)
calculate TET volume
Definition: fem_tools.c:298
FTensor::Index< 'n', SPACE_DIM > n
FTensor::Index< 'm', SPACE_DIM > m
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
FTensor::Tensor1< double, SPACE_DIM > normal(FTensor::Tensor1< T1, 3 > &t_coords, FTensor::Tensor1< T2, SPACE_DIM > &t_disp)
Definition: ContactOps.hpp:65
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
static std::map< long int, MatrixDouble > mapRefCoords
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
VectorBoundedArray< double, 3 > VectorDouble3
Definition: Types.hpp:92
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:40
implementation of Data Operators for Forces and Sources
Definition: MoFEM.hpp:24
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:1210
CoreTmp< 0 > Core
Definition: Core.hpp:1086
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:135
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
Definition: Templates.hpp:1199
double w(const double g, const double t)
const double r
rate factor
double rho
Definition: plastic.cpp:98
constexpr auto field_name
#define QUAD_3D_TABLE_SIZE
Definition: quad.h:186
static QUAD *const QUAD_3D_TABLE[]
Definition: quad.h:187
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)
MoFEMErrorCode iNtegrate(DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
boost::shared_ptr< HookeElement::DataAtIntegrationPts > dataAtPts
MoFEMErrorCode aSsemble(DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
MoFEMErrorCode iNtegrate(DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
boost::shared_ptr< HookeElement::DataAtIntegrationPts > dataAtPts
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)
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
MoFEMErrorCode aSsemble(DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
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)
boost::shared_ptr< MatrixDouble > hG
spatial deformation gradient
boost::ptr_vector< NeumannForcesSurface::MethodForAnalyticalForce > & analyticalForceOp
boost::ptr_vector< MethodForForceScaling > & methodsOp
VectorDouble sNrm
Length of the normal vector.
boost::shared_ptr< MatrixDouble > HG
spatial deformation gradient
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
CrackFrontSingularBase< FirendVolumeOnSide, VolumeElementForcesAndSourcesCore > volSideFe
boost::ptr_vector< MethodForForceScaling > & methodsOp
CrackFrontSingularBase< VolumeElementForcesAndSourcesCoreOnSide, VolumeElementForcesAndSourcesCore > volSideFe
VectorDouble sNrm
Length of the normal vector.
boost::ptr_vector< NeumannForcesSurface::MethodForAnalyticalForce > & analyticalForceOp
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)
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
MoFEMErrorCode doWork(int row_side, EntityType row_type, HookeElement::EntData &row_data)
boost::shared_ptr< MatrixDouble > matGradPosAtPtsPtr
boost::shared_ptr< CrackFrontElement > feSingularPtr
boost::shared_ptr< VectorDouble > rhoAtGaussPtsPtr
boost::shared_ptr< MatrixDouble > rhoGradAtGaussPtsPtr
boost::shared_ptr< MatrixDouble > matCoordsPtr
boost::shared_ptr< MatrixDouble > rhoGradGradAtGaussPtsPtr
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)
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
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)
moab::Interface & postProcMesh
MatrixDouble & singularDisp
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
bool & singularElement
std::vector< EntityHandle > & mapGaussPts
boost::shared_ptr< MatrixDouble > H
OpRhsBoneExplicitDerivariveWithHooke(HookeElement::DataAtIntegrationPts &common_data, boost::shared_ptr< VectorDouble > rho_at_gauss_pts, boost::shared_ptr< MatrixDouble > rho_grad_at_gauss_pts, Range tets_in_block, const double rho_0, const double bone_n, Range *forces_only_on_entities_row=NULL)
MoFEMErrorCode doWork(int row_side, EntityType row_type, DataForcesAndSourcesCore::EntData &row_data)
MoFEMErrorCode doWork(int row_side, EntityType row_type, DataForcesAndSourcesCore::EntData &row_data)
std::vector< EntityHandle > & mapGaussPts
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
static MoFEMErrorCode applyScale(const FEMethod *fe, boost::ptr_vector< MethodForForceScaling > &methods_op, VectorDouble &nf)
Managing BitRefLevels.
Deprecated interface functions.
EntityHandle getFEEntityHandle() const
Return finite element entity handle.
boost::shared_ptr< const NumeredEntFiniteElement > getNumeredEntFiniteElementPtr() const
Return raw pointer to NumeredEntFiniteElement.
const FEMethod * getFEMethod() const
Return raw pointer to Finite Element Method object.
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
Mesh refinement interface.
MoFEMErrorCode refineTets(const EntityHandle meshset, const BitRefLevel &bit, int verb=QUIET, const bool debug=false)
refine TET in the meshset
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
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
operator calculating deformation gradient
int order
Definition: quad.h:28
int npoints
Definition: quad.h:29