v0.15.5
Loading...
Searching...
No Matches
EshelbianTopologicalDerivativeOperators.cpp
Go to the documentation of this file.
1/**
2 * @file EshelbianTopologicalDerivative.cpp
3 * @brief
4 * @version 0.1
5 * @date 2026-02-11
6 *
7 * @copyright Copyright (c) 2026
8 *
9 */
10
11// #define SINGULARITY
12// #include <MoFEM.hpp>
13// using namespace MoFEM;
14
15// #include <EshelbianPlasticity.hpp>
16// #include <EshelbianAux.hpp>
17// #include <TSElasticPostStep.hpp>
18
19// #include <boost/math/constants/constants.hpp>
20// #include <boost/math/special_functions/lambert_w.hpp>
21
22// #include <EshelbianTopologicalDerivative.hpp>
23// #include <ObjectiveFunctionData.hpp>
24
25using namespace EshelbianPlasticity;
26using namespace ShapeOptimization;
27
28namespace EshelbianPlasticity {
29
30// dr/dX topological derivative for Eshelbian plasticity model
31
32template <typename AssembleOp>
33struct OpAssembleTopologicalObjectiveDerivativeImplBase : public AssembleOp {
34 using OP = AssembleOp;
35
36 OpAssembleTopologicalObjectiveDerivativeImplBase(
37 const std::string &field_name,
38 boost::shared_ptr<DataAtIntegrationPts> data_ptr,
39 boost::shared_ptr<TopologicalData> topo_ptr,
40 boost::shared_ptr<double> J_ptr, SmartPetscObj<Vec> assemble_vec,
41 Tag topo_tag)
42 : OP(field_name, data_ptr, OP::OPROW), JPtr(J_ptr),
43 assembleVec(assemble_vec), topoTag(topo_tag), topoData(topo_ptr) {}
44
45 MoFEMErrorCode assemble(EntData &data) override {
46 MoFEMFunctionBegin;
47 if (assembleVec) {
48 double *vec_ptr = OP::nF.data().data();
49 const int nb_dofs = data.getIndices().size();
50 int *ind_ptr = data.getIndices().data().data();
51 CHKERR VecSetValues(assembleVec, nb_dofs, ind_ptr, vec_ptr, ADD_VALUES);
52 }
53 if (locJ) {
54 *JPtr += locJ;
55 }
56 if (topoTag) {
57 const auto field_ents = data.getFieldEntities();
58 std::vector<EntityHandle> ents(field_ents.size());
59 std::transform(field_ents.begin(), field_ents.end(), ents.begin(),
60 [](const auto *fe) { return fe->getEnt(); });
61 if (field_ents.empty())
62 MoFEMFunctionReturnHot(0);
63 if (type_from_handle(ents[0]) != MBVERTEX)
64 MoFEMFunctionReturnHot(0);
65 auto &moab = OP::getMoab();
66 VectorDouble topo_values(OP::nF.size());
67 CHKERR moab.tag_set_data(topoTag, ents.data(), ents.size(),
68 topo_values.data().data());
69 noalias(topo_values) += OP::nF;
70 CHKERR moab.tag_set_data(topoTag, ents.data(), ents.size(),
71 OP::nF.data().data());
72 }
73 MoFEMFunctionReturn(0);
74 }
75
76protected:
77 double locJ;
78 boost::shared_ptr<double> JPtr;
79 SmartPetscObj<Vec> assembleVec;
80 Tag topoTag;
81 boost::shared_ptr<TopologicalData> topoData;
82};
83
90
96
97struct OpAssembleBrokenFaceTopologicalDerivativeImplBase
98 : public FormsIntegrators<FaceUserDataOperator>::Assembly<A>::OpBrokenBase {
99
100 using OP = typename FormsIntegrators<FaceUserDataOperator>::Assembly<
101 A>::OpBrokenBase;
102
103 OpAssembleBrokenFaceTopologicalDerivativeImplBase(
104 boost::shared_ptr<std::vector<BrokenBaseSideData>> broken_base_side_data,
105 boost::shared_ptr<TopologicalData> topo_ptr,
106 boost::shared_ptr<double> J_ptr,
107 SmartPetscObj<Vec> assemble_vec,
108 Tag topo_tag,
109 boost::shared_ptr<Range> ents_ptr = nullptr)
110 : OP(broken_base_side_data, ents_ptr), JPtr(J_ptr),
111 assembleVec(assemble_vec), topoTag(topo_tag),
112 topoData(topo_ptr) {}
113
114 MoFEMErrorCode aSsemble(EntData &data) override {
115 MoFEMFunctionBegin;
116 if (assembleVec) {
117 double *vec_ptr = OP::locF.data().data();
118 const int nb_dofs = data.getIndices().size();
119 int *ind_ptr = data.getIndices().data().data();
120 CHKERR VecSetValues(assembleVec, nb_dofs, ind_ptr, vec_ptr, ADD_VALUES);
121 }
122 if (topoTag) {
123 const auto field_ents = data.getFieldEntities();
124 std::vector<EntityHandle> ents(field_ents.size());
125 std::transform(field_ents.begin(), field_ents.end(), ents.begin(),
126 [](const auto *fe) { return fe->getEnt(); });
127 if (field_ents.empty())
128 MoFEMFunctionReturnHot(0);
129 if (type_from_handle(ents[0]) != MBVERTEX)
130 MoFEMFunctionReturnHot(0);
131 auto &moab = getMoab();
132 VectorDouble topo_values(OP::locF.size());
133 CHKERR moab.tag_set_data(topoTag, ents.data(), ents.size(),
134 topo_values.data().data());
135 topo_values += OP::locF;
136 CHKERR moab.tag_set_data(topoTag, ents.data(), ents.size(),
137 OP::locF.data().data());
138 }
139 MoFEMFunctionReturn(0);
140 }
141
142protected:
143 boost::shared_ptr<double> JPtr;
144 SmartPetscObj<Vec> assembleVec;
145 Tag topoTag;
146 boost::shared_ptr<TopologicalData> topoData;
147};
148
149struct OpJ_dPImpl : public OpAssembleVolumeTopologicalDerivativeImpl {
150 using OpAssembleVolumeTopologicalDerivativeImpl::
151 OpAssembleVolumeTopologicalDerivativeImpl;
152
153 MoFEMErrorCode integrate(EntData &data);
154};
155
156struct OpJ_dBubbleImpl : public OpAssembleVolumeTopologicalDerivativeImpl {
157 using OpAssembleVolumeTopologicalDerivativeImpl::
158 OpAssembleVolumeTopologicalDerivativeImpl;
159
160 MoFEMErrorCode integrate(EntData &data);
161};
162
163struct OpJ_dwImpl : public OpAssembleVolumeTopologicalDerivativeImpl {
164 using OpAssembleVolumeTopologicalDerivativeImpl::
165 OpAssembleVolumeTopologicalDerivativeImpl;
166
167 MoFEMErrorCode integrate(EntData &data);
168};
169
170struct dJ_duGammaImpl : public OpAssembleFaceTopologicalDerivativeImpl {
171 using OpAssembleFaceTopologicalDerivativeImpl::
172 OpAssembleFaceTopologicalDerivativeImpl;
173
174 MoFEMErrorCode integrate(EntData &data);
175};
176
177struct dJ_dTractionImpl
178 : public OpAssembleBrokenFaceTopologicalDerivativeImplBase {
179 using OpAssembleBrokenFaceTopologicalDerivativeImplBase::
180 OpAssembleBrokenFaceTopologicalDerivativeImplBase;
181
182 MoFEMErrorCode iNtegrate(EntData &data);
183};
184
185struct OpJ_dP_no_streachImpl
186 : public OpAssembleVolumeTopologicalDerivativeImpl {
187 using OpAssembleVolumeTopologicalDerivativeImpl::
188 OpAssembleVolumeTopologicalDerivativeImpl;
189
190 MoFEMErrorCode integrate(EntData &data);
191};
192
193struct OpJ_dBubble_no_streachImpl
194 : public OpAssembleVolumeTopologicalDerivativeImpl {
195 using OpAssembleVolumeTopologicalDerivativeImpl::
196 OpAssembleVolumeTopologicalDerivativeImpl;
197
198 MoFEMErrorCode integrate(EntData &data);
199};
200
201struct OpJ_dUImpl : public OpAssembleVolumeTopologicalDerivativeImpl {
202 using OpAssembleVolumeTopologicalDerivativeImpl::
203 OpAssembleVolumeTopologicalDerivativeImpl;
204
205 MoFEMErrorCode integrate(EntData &data);
206};
207
208struct OpInteriorJImpl : public OpAssembleVolumeTopologicalDerivativeImpl {
209 using OpAssembleVolumeTopologicalDerivativeImpl::
210 OpAssembleVolumeTopologicalDerivativeImpl;
211
212 MoFEMErrorCode integrate(EntData &data);
213};
214
215struct OpTopologicalObjectivePythonImpl
216 : public ForcesAndSourcesCore::UserDataOperator {
217 using OP = ForcesAndSourcesCore::UserDataOperator;
218
219 OpTopologicalObjectivePythonImpl(
220 boost::shared_ptr<DataAtIntegrationPts> data_at_pts_ptr,
221 boost::shared_ptr<TopologicalData> topo_p,
222 boost::shared_ptr<ObjectiveFunctionData> python_ptr)
223 : OP(NOSPACE, OP::OPSPACE), dataAtPts(data_at_pts_ptr), topoData(topo_p),
224 pythonPtr(python_ptr) {}
225
226 MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
227
228private:
229 boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
230 boost::shared_ptr<TopologicalData> topoData;
231 boost::shared_ptr<ObjectiveFunctionData> pythonPtr;
232};
233
234MoFEMErrorCode OpTopologicalObjectivePythonImpl::doWork(int side,
235 EntityType type,
236 EntData &data) {
237 MoFEMFunctionBegin;
238
239#ifndef NDEBUG
240 if (!dataAtPts)
241 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
242 "DataAtIntegrationPts pointer is null");
243 if (!topoData)
244 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
245 "Topological data pointer is null");
246 if (!pythonPtr)
247 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
248 "ObjectiveFunctionData pointer is null");
249#endif // NDEBUG
250
251 constexpr int size_full = SPACE_DIM * SPACE_DIM;
252 const int nb_gauss_pts = getGaussPts().size2();
253 if (!nb_gauss_pts)
254 MoFEMFunctionReturnHot(0);
255
256 auto stress_full_ptr = boost::make_shared<MatrixDouble>();
257 stress_full_ptr->resize(size_full, nb_gauss_pts, false);
258 stress_full_ptr->clear();
259 auto strain_full_ptr = boost::make_shared<MatrixDouble>();
260 strain_full_ptr->resize(size_full, nb_gauss_pts, false);
261 strain_full_ptr->clear();
262
263 auto t_stress = getFTensor2FromMat<3, 3>(*stress_full_ptr);
264 auto t_strain = getFTensor2FromMat<3, 3>(*strain_full_ptr);
265
266 auto t_P = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(dataAtPts->approxPAtPts);
267 auto t_log_u =
268 getFTensor2SymmetricFromMat<3>(dataAtPts->logStretchTensorAtPts);
269
270 auto next = [&]() {
271 ++t_stress;
272 ++t_strain;
273 ++t_P;
274 ++t_log_u;
275 };
276
277 FTENSOR_INDEX(SPACE_DIM, i);
278 FTENSOR_INDEX(SPACE_DIM, j);
279
280 for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
281 // we have to handle all variants, that will render how the Jacobian
282 // gradient is evaluated.
283 t_stress(i, j) = t_P(i, j);
284 t_strain(i, j) = t_log_u(i, j);
285 next();
286 }
287
288 auto &coords = OP::getCoordsAtGaussPts();
289 CHKERR pythonPtr->evalInteriorObjectiveFunction(
290 coords, dataAtPts->getSmallWL2AtPts(), stress_full_ptr, strain_full_ptr,
291 topoData->getObjAtPts(), false);
292
293 CHKERR pythonPtr->evalInteriorObjectiveGradientStrain(
294 coords, dataAtPts->getSmallWL2AtPts(), stress_full_ptr, strain_full_ptr,
295 topoData->getObjDStrainAtPts(), false);
296
297 CHKERR pythonPtr->evalInteriorObjectiveGradientU(
298 coords, dataAtPts->getSmallWL2AtPts(), stress_full_ptr, strain_full_ptr,
299 topoData->getObjDDisplacementAtPts(), false);
300
301 CHKERR pythonPtr->evalInteriorObjectiveGradientStress(
302 coords, dataAtPts->getSmallWL2AtPts(), stress_full_ptr, strain_full_ptr,
303 topoData->getObjDStressAtPts(), false);
304
305 MoFEMFunctionReturn(0);
306}
307
308MoFEMErrorCode OpInteriorJImpl::integrate(EntData &data) {
309 MoFEMFunctionBegin;
310
311#ifndef NDEBUG
312 if (!topoData)
313 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
314 "Topological data pointer is null");
315#endif // NDEBUG
316
317 const int nb_dofs = data.getIndices().size();
318 if (!nb_dofs)
319 MoFEMFunctionReturnHot(0);
320
321 const int nb_integration_pts = data.getN().size1();
322 auto obj_at_pts = topoData->getObjAtPts();
323
324#ifndef NDEBUG
325 if (obj_at_pts->size() != static_cast<size_t>(nb_integration_pts))
326 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
327 "objAtPts size (%d) != nb integration points (%d)",
328 static_cast<int>(obj_at_pts->size()), nb_integration_pts);
329 if (topoData->detJacobianAtPts.size() !=
330 static_cast<size_t>(nb_integration_pts))
331 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
332 "detJacobianAtPts size (%d) != nb integration points (%d)",
333 static_cast<int>(topoData->detJacobianAtPts.size()),
334 nb_integration_pts);
335 if (topoData->invJacobianAtPtr.size2() !=
336 static_cast<size_t>(nb_integration_pts))
337 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
338 "invJacobianAtPtr columns (%d) != nb integration points (%d)",
339 static_cast<int>(topoData->invJacobianAtPtr.size2()),
340 nb_integration_pts);
341 if (topoData->invJacobianAtPtr.size1() != SPACE_DIM * SPACE_DIM)
342 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
343 "invJacobianAtPtr row size is %d, expected %d",
344 static_cast<int>(topoData->invJacobianAtPtr.size1()),
345 SPACE_DIM * SPACE_DIM);
346#endif // NDEBUG
347
348 const auto v = getVolume();
349 auto t_w = getFTensor0IntegrationWeight();
350 auto t_obj = getFTensor0FromVec(*obj_at_pts);
351 auto t_det = getFTensor0FromVec(topoData->detJacobianAtPts);
352 auto t_inv_jac =
353 getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(topoData->invJacobianAtPtr);
354
355 const int nb_base_functions = data.getN().size2();
356 auto t_base_diff = data.getFTensor1DiffN<3>();
357
358 FTENSOR_INDEX(SPACE_DIM, i);
359 FTENSOR_INDEX(SPACE_DIM, j);
360
361 auto get_ftensor1 = [](auto &v) {
362 return FTensor::Tensor1<FTensor::PackPtr<double *, SPACE_DIM>, SPACE_DIM>(
363 &v[0], &v[1], &v[2]);
364 };
365
366 locJ = 0;
367
368 for (int gg = 0; gg != nb_integration_pts; ++gg) {
369 const double a = v * t_w;
370
371 locJ += a * t_obj;
372
373 FTensor::Tensor2<double, SPACE_DIM, SPACE_DIM> t_cof;
374 t_cof(i, j) = t_det * t_inv_jac(j, i);
375
376 auto t_nf = get_ftensor1(nF);
377 int bb = 0;
378 for (; bb != nb_dofs / SPACE_DIM; ++bb) {
379 t_nf(i) += a * t_obj * t_cof(i, j) * t_base_diff(j);
380 ++t_nf;
381 ++t_base_diff;
382 }
383 for (; bb != nb_base_functions; ++bb)
384 ++t_base_diff;
385
386 ++t_w;
387 ++t_obj;
388 ++t_det;
389 ++t_inv_jac;
390 }
391
392 MoFEMFunctionReturn(0);
393}
394
395MoFEMErrorCode OpJ_dPImpl::integrate(EntData &data) {
396 MoFEMFunctionBegin;
397
398#ifndef NDEBUG
399 if (!topoData)
400 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
401 "Topological data pointer is null");
402#endif // NDEBUG
403
404 const int nb_dofs = data.getIndices().size();
405 if (!nb_dofs)
406 MoFEMFunctionReturnHot(0);
407
408 const int nb_integration_pts = data.getN().size1();
409
410#ifndef NDEBUG
411 const int stress_rows = topoData->objDStressAtPts.size1();
412 const int stress_cols = topoData->objDStressAtPts.size2();
413 if (stress_cols != nb_integration_pts)
414 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
415 "objDStressAtPts columns (%d) != nb integration points (%d)",
416 stress_cols, nb_integration_pts);
417
418 if (stress_rows != SPACE_DIM * SPACE_DIM) {
419 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
420 "objDStressAtPts row size is %d, expected %d (3x3 Piola gradient)",
421 stress_rows, SPACE_DIM * SPACE_DIM);
422 }
423#endif // NDEBUG
424
425 const auto v = getVolume();
426 auto t_w = getFTensor0IntegrationWeight();
427 const int nb_base_functions = data.getN().size2() / SPACE_DIM;
428 auto t_row_base_fun = data.getFTensor1N<SPACE_DIM>();
429
430 FTENSOR_INDEX(SPACE_DIM, i);
431 FTENSOR_INDEX(SPACE_DIM, j);
432
433 auto get_ftensor1 = [](auto &v) {
434 return FTensor::Tensor1<FTensor::PackPtr<double *, SPACE_DIM>, SPACE_DIM>(
435 &v[0], &v[1], &v[2]);
436 };
437
438 auto t_obj_dP =
439 getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(topoData->objDStressAtPts);
440
441 for (int gg = 0; gg != nb_integration_pts; ++gg) {
442 const double a = v * t_w;
443 auto t_nf = get_ftensor1(nF);
444
445 int bb = 0;
446 for (; bb != nb_dofs / SPACE_DIM; ++bb) {
447 t_nf(i) += a * t_row_base_fun(j) * t_obj_dP(i, j);
448 ++t_nf;
449 ++t_row_base_fun;
450 }
451 for (; bb != nb_base_functions; ++bb)
452 ++t_row_base_fun;
453
454 ++t_w;
455 ++t_obj_dP;
456 }
457
458 MoFEMFunctionReturn(0);
459}
460
461MoFEMErrorCode OpJ_dBubbleImpl::integrate(EntData &data) {
462 MoFEMFunctionBegin;
463
464#ifndef NDEBUG
465 if (!topoData)
466 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
467 "Topological data pointer is null");
468#endif // NDEBUG
469
470 const int nb_dofs = data.getIndices().size();
471 if (!nb_dofs)
472 MoFEMFunctionReturnHot(0);
473
474 const int nb_integration_pts = data.getN().size1();
475
476#ifndef NDEBUG
477 const int stress_rows = topoData->objDStressAtPts.size1();
478 const int stress_cols = topoData->objDStressAtPts.size2();
479 if (stress_cols != nb_integration_pts)
480 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
481 "objDStressAtPts columns (%d) != nb integration points (%d)",
482 stress_cols, nb_integration_pts);
483
484 if (stress_rows != SPACE_DIM * SPACE_DIM) {
485 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
486 "objDStressAtPts row size is %d, expected %d (3x3 Piola gradient)",
487 stress_rows, SPACE_DIM * SPACE_DIM);
488 }
489#endif // NDEBUG
490
491 const auto v = getVolume();
492 auto t_w = getFTensor0IntegrationWeight();
493 const int nb_base_functions = data.getN().size2() / (SPACE_DIM * SPACE_DIM);
494 auto t_row_base_fun = data.getFTensor2N<SPACE_DIM, SPACE_DIM>();
495 auto t_obj_dP =
496 getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(topoData->objDStressAtPts);
497
498 FTENSOR_INDEX(SPACE_DIM, i);
499 FTENSOR_INDEX(SPACE_DIM, j);
500
501 auto get_ftensor0 = [](auto &v) {
502 return FTensor::Tensor0<FTensor::PackPtr<double *, 1>>(&v[0]);
503 };
504
505 for (int gg = 0; gg != nb_integration_pts; ++gg) {
506 const double a = v * t_w;
507 auto t_nf = get_ftensor0(nF);
508
509 int bb = 0;
510 for (; bb != nb_dofs; ++bb) {
511 t_nf += a * t_row_base_fun(i, j) * t_obj_dP(i, j);
512 ++t_nf;
513 ++t_row_base_fun;
514 }
515 for (; bb != nb_base_functions; ++bb)
516 ++t_row_base_fun;
517
518 ++t_w;
519 ++t_obj_dP;
520 }
521
522 MoFEMFunctionReturn(0);
523}
524
525MoFEMErrorCode OpJ_dwImpl::integrate(EntData &data) {
526 MoFEMFunctionBegin;
527
528#ifndef NDEBUG
529 if (!topoData)
530 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
531 "Topological data pointer is null");
532#endif // NDEBUG
533
534 const int nb_dofs = data.getIndices().size();
535 if (!nb_dofs)
536 MoFEMFunctionReturnHot(0);
537
538 const int nb_integration_pts = data.getN().size1();
539
540#ifndef NDEBUG
541 const int disp_rows = topoData->objDDisplacementAtPts.size1();
542 const int disp_cols = topoData->objDDisplacementAtPts.size2();
543 if (disp_cols != nb_integration_pts)
544 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
545 "objDDisplacementAtPts columns (%d) != nb integration points (%d)",
546 disp_cols, nb_integration_pts);
547
548 if (disp_rows != SPACE_DIM) {
549 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
550 "objDDisplacementAtPts row size is %d, expected %d", disp_rows,
551 SPACE_DIM);
552 }
553#endif // NDEBUG
554
555 const auto v = getVolume();
556 auto t_w = getFTensor0IntegrationWeight();
557 const int nb_base_functions = data.getN().size2();
558 auto t_row_base_fun = data.getFTensor0N();
559 auto t_obj_dw =
560 getFTensor1FromMat<SPACE_DIM>(topoData->objDDisplacementAtPts);
561
562 FTENSOR_INDEX(SPACE_DIM, i);
563
564 auto get_ftensor1 = [](auto &v) {
565 return FTensor::Tensor1<FTensor::PackPtr<double *, SPACE_DIM>, SPACE_DIM>(
566 &v[0], &v[1], &v[2]);
567 };
568
569 for (int gg = 0; gg != nb_integration_pts; ++gg) {
570 const double a = v * t_w;
571 auto t_nf = get_ftensor1(nF);
572
573 int bb = 0;
574 for (; bb != nb_dofs / SPACE_DIM; ++bb) {
575 t_nf(i) += a * t_row_base_fun * t_obj_dw(i);
576 ++t_nf;
577 ++t_row_base_fun;
578 }
579 for (; bb != nb_base_functions; ++bb)
580 ++t_row_base_fun;
581
582 ++t_w;
583 ++t_obj_dw;
584 }
585
586 MoFEMFunctionReturn(0);
587}
588
589MoFEMErrorCode dJ_duGammaImpl::integrate(EntData &data) {
590 MoFEMFunctionBegin;
591
592#ifndef NDEBUG
593 if (!topoData)
594 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
595 "Topological data pointer is null");
596#endif // NDEBUG
597
598 const int nb_dofs = data.getIndices().size();
599 if (!nb_dofs)
600 MoFEMFunctionReturnHot(0);
601
602 const int nb_integration_pts = data.getN().size1();
603
604#ifndef NDEBUG
605 const int disp_rows = topoData->objDDisplacementAtPts.size1();
606 const int disp_cols = topoData->objDDisplacementAtPts.size2();
607 if (disp_cols != nb_integration_pts)
608 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
609 "objDDisplacementAtPts columns (%d) != nb integration points (%d)",
610 disp_cols, nb_integration_pts);
611
612 if (disp_rows != SPACE_DIM) {
613 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
614 "objDDisplacementAtPts row size is %d, expected %d", disp_rows,
615 SPACE_DIM);
616 }
617#endif // NDEBUG
618
619 auto t_w = getFTensor0IntegrationWeight();
620 const int nb_base_functions = data.getN().size2();
621 auto t_row_base_fun = data.getFTensor0N();
622 auto t_obj_du_gamma =
623 getFTensor1FromMat<SPACE_DIM>(topoData->objDDisplacementAtPts);
624
625 FTENSOR_INDEX(SPACE_DIM, i);
626
627 auto get_ftensor1 = [](auto &v) {
628 return FTensor::Tensor1<FTensor::PackPtr<double *, SPACE_DIM>, SPACE_DIM>(
629 &v[0], &v[1], &v[2]);
630 };
631
632 for (int gg = 0; gg != nb_integration_pts; ++gg) {
633 const double a = t_w * getMeasure();
634 auto t_nf = get_ftensor1(nF);
635
636 int bb = 0;
637 for (; bb != nb_dofs / SPACE_DIM; ++bb) {
638 t_nf(i) += a * t_row_base_fun * t_obj_du_gamma(i);
639 ++t_nf;
640 ++t_row_base_fun;
641 }
642 for (; bb != nb_base_functions; ++bb)
643 ++t_row_base_fun;
644
645 ++t_w;
646 ++t_obj_du_gamma;
647 }
648
649 MoFEMFunctionReturn(0);
650}
651
652MoFEMErrorCode dJ_dTractionImpl::iNtegrate(EntData &data) {
653 MoFEMFunctionBegin;
654
655#ifndef NDEBUG
656 if (!topoData)
657 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
658 "Topological data pointer is null");
659#endif // NDEBUG
660
661 const int nb_dofs = data.getIndices().size();
662 if (!nb_dofs)
663 MoFEMFunctionReturnHot(0);
664
665 const int nb_integration_pts = OP::getGaussPts().size2();
666
667#ifndef NDEBUG
668 const int traction_rows = topoData->objDDisplacementAtPts.size1();
669 const int traction_cols = topoData->objDDisplacementAtPts.size2();
670 if (traction_cols != nb_integration_pts)
671 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
672 "objDDisplacementAtPts columns (%d) != nb integration points (%d)",
673 traction_cols, nb_integration_pts);
674
675 if (traction_rows != SPACE_DIM) {
676 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
677 "objDDisplacementAtPts row size is %d, expected %d", traction_rows,
678 SPACE_DIM);
679 }
680#endif // NDEBUG
681
682 auto t_normal = OP::getFTensor1NormalsAtGaussPts();
683 auto t_w = OP::getFTensor0IntegrationWeight();
684 const int nb_base_functions = data.getN().size2() / SPACE_DIM;
685 auto t_row_base_fun = data.getFTensor1N<SPACE_DIM>();
686 auto t_obj_dtraction =
687 getFTensor1FromMat<SPACE_DIM>(topoData->objDDisplacementAtPts);
688
689 FTENSOR_INDEX(SPACE_DIM, i);
690 FTENSOR_INDEX(SPACE_DIM, j);
691
692 for (int gg = 0; gg != nb_integration_pts; ++gg) {
693 auto t_nf = getFTensor1FromPtr<SPACE_DIM>(&*OP::locF.begin());
694 int bb = 0;
695 for (; bb != nb_dofs / SPACE_DIM; ++bb) {
696 t_nf(i) +=
697 t_w * (t_row_base_fun(j) * t_normal(j)) * t_obj_dtraction(i) * 0.5;
698 ++t_nf;
699 ++t_row_base_fun;
700 }
701 for (; bb != nb_base_functions; ++bb)
702 ++t_row_base_fun;
703
704 ++t_w;
705 ++t_normal;
706 ++t_obj_dtraction;
707 }
708
709 MoFEMFunctionReturn(0);
710}
711
712MoFEMErrorCode OpJ_dP_no_streachImpl::integrate(EntData &data) {
713 MoFEMFunctionBegin;
714
715#ifndef NDEBUG
716 if (!topoData)
717 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
718 "Topological data pointer is null");
719 if (!dataAtPts)
720 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
721 "DataAtIntegrationPts pointer is null");
722#endif // NDEBUG
723
724 const int nb_dofs = data.getIndices().size();
725 if (!nb_dofs)
726 MoFEMFunctionReturnHot(0);
727
728 const int nb_integration_pts = data.getN().size1();
729
730#ifndef NDEBUG
731 const int strain_rows = topoData->objDStrainAtPts.size1();
732 const int strain_cols = topoData->objDStrainAtPts.size2();
733 if (strain_cols != nb_integration_pts)
734 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
735 "objDStrainAtPts columns (%d) != nb integration points (%d)",
736 strain_cols, nb_integration_pts);
737
738 if (strain_rows != size_symm) {
739 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
740 "objDStrainAtPts row size is %d, expected %d", strain_rows,
741 size_symm);
742 }
743#endif // NDEBUG
744
745 auto &mat_d = dataAtPts->matD;
746#ifndef NDEBUG
747 if (mat_d.size1() != size_symm || mat_d.size2() != size_symm) {
748 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
749 "wrong matD size, should be %d by %d but is %d by %d", size_symm,
750 size_symm, (int)mat_d.size1(), (int)mat_d.size2());
751 }
752#endif // NDEBUG
753
754 const auto v = getVolume();
755 auto t_w = getFTensor0IntegrationWeight();
756 const int nb_base_functions = data.getN().size2() / SPACE_DIM;
757 auto t_row_base_fun = data.getFTensor1N<SPACE_DIM>();
758
759 auto t_obj_dstrain =
760 getFTensor2SymmetricFromMat<SPACE_DIM>(topoData->objDStrainAtPts);
761 auto t_D =
762 getFTensor4DdgFromPtr<SPACE_DIM, SPACE_DIM, 0>(mat_d.data().data());
763
764 FTENSOR_INDEX(SPACE_DIM, i);
765 FTENSOR_INDEX(SPACE_DIM, j);
766 FTENSOR_INDEX(SPACE_DIM, k);
767 FTENSOR_INDEX(SPACE_DIM, l);
768
769 auto get_ftensor1 = [](auto &v) {
770 return FTensor::Tensor1<FTensor::PackPtr<double *, SPACE_DIM>, SPACE_DIM>(
771 &v[0], &v[1], &v[2]);
772 };
773
774 for (int gg = 0; gg != nb_integration_pts; ++gg) {
775 const double a = v * t_w;
776 auto t_nf = get_ftensor1(nF);
777
778 FTensor::Tensor2_symmetric<double, SPACE_DIM> t_dj_dp_no_streach;
779 t_dj_dp_no_streach(i, j) = t_obj_dstrain(k, l) * t_D(k, l, i, j);
780
781 int bb = 0;
782 for (; bb != nb_dofs / SPACE_DIM; ++bb) {
783 t_nf(i) += a * t_row_base_fun(j) * t_dj_dp_no_streach(i, j);
784 ++t_nf;
785 ++t_row_base_fun;
786 }
787 for (; bb != nb_base_functions; ++bb)
788 ++t_row_base_fun;
789
790 ++t_w;
791 ++t_obj_dstrain;
792 }
793
794 MoFEMFunctionReturn(0);
795}
796
797MoFEMErrorCode OpJ_dBubble_no_streachImpl::integrate(EntData &data) {
798 MoFEMFunctionBegin;
799
800#ifndef NDEBUG
801 if (!topoData)
802 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
803 "Topological data pointer is null");
804 if (!dataAtPts)
805 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
806 "DataAtIntegrationPts pointer is null");
807#endif // NDEBUG
808
809 const int nb_dofs = data.getIndices().size();
810 if (!nb_dofs)
811 MoFEMFunctionReturnHot(0);
812
813 const int nb_integration_pts = data.getN().size1();
814
815#ifndef NDEBUG
816 const int strain_rows = topoData->objDStrainAtPts.size1();
817 const int strain_cols = topoData->objDStrainAtPts.size2();
818 if (strain_cols != nb_integration_pts)
819 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
820 "objDStrainAtPts columns (%d) != nb integration points (%d)",
821 strain_cols, nb_integration_pts);
822
823 if (strain_rows != size_symm) {
824 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
825 "objDStrainAtPts row size is %d, expected %d", strain_rows,
826 size_symm);
827 }
828#endif // NDEBUG
829
830 auto &mat_d = dataAtPts->matD;
831#ifndef NDEBUG
832 if (mat_d.size1() != size_symm || mat_d.size2() != size_symm) {
833 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
834 "wrong matD size, should be %d by %d but is %d by %d", size_symm,
835 size_symm, (int)mat_d.size1(), (int)mat_d.size2());
836 }
837#endif // NDEBUG
838
839 const auto v = getVolume();
840 auto t_w = getFTensor0IntegrationWeight();
841 const int nb_base_functions = data.getN().size2() / (SPACE_DIM * SPACE_DIM);
842 auto t_row_base_fun = data.getFTensor2N<SPACE_DIM, SPACE_DIM>();
843
844 auto t_obj_dstrain =
845 getFTensor2SymmetricFromMat<SPACE_DIM>(topoData->objDStrainAtPts);
846 auto t_D =
847 getFTensor4DdgFromPtr<SPACE_DIM, SPACE_DIM, 0>(mat_d.data().data());
848
849 FTENSOR_INDEX(SPACE_DIM, i);
850 FTENSOR_INDEX(SPACE_DIM, j);
851 FTENSOR_INDEX(SPACE_DIM, k);
852 FTENSOR_INDEX(SPACE_DIM, l);
853
854 auto get_ftensor0 = [](auto &v) {
855 return FTensor::Tensor0<FTensor::PackPtr<double *, 1>>(&v[0]);
856 };
857
858 for (int gg = 0; gg != nb_integration_pts; ++gg) {
859 const double a = v * t_w;
860 auto t_nf = get_ftensor0(nF);
861
862 FTensor::Tensor2_symmetric<double, SPACE_DIM> t_dj_dp_no_streach;
863 t_dj_dp_no_streach(i, j) = t_obj_dstrain(k, l) * t_D(k, l, i, j);
864
865 int bb = 0;
866 for (; bb != nb_dofs; ++bb) {
867 t_nf += a * t_row_base_fun(i, j) * t_dj_dp_no_streach(i, j);
868 ++t_nf;
869 ++t_row_base_fun;
870 }
871 for (; bb != nb_base_functions; ++bb)
872 ++t_row_base_fun;
873
874 ++t_w;
875 ++t_obj_dstrain;
876 }
877
878 MoFEMFunctionReturn(0);
879}
880
881MoFEMErrorCode OpJ_dUImpl::integrate(EntData &data) {
882 MoFEMFunctionBegin;
883
884#ifndef NDEBUG
885 if (!topoData)
886 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
887 "Topological data pointer is null");
888#endif // NDEBUG
889
890 const int nb_dofs = data.getIndices().size();
891 if (!nb_dofs)
892 MoFEMFunctionReturnHot(0);
893
894 const int nb_integration_pts = data.getN().size1();
895
896#ifndef NDEBUG
897 const int strain_rows = topoData->objDStrainAtPts.size1();
898 const int strain_cols = topoData->objDStrainAtPts.size2();
899 if (strain_cols != nb_integration_pts)
900 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
901 "objDStrainAtPts columns (%d) != nb integration points (%d)",
902 strain_cols, nb_integration_pts);
903
904 if (strain_rows != size_symm) {
905 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
906 "objDStrainAtPts row size is %d, expected %d", strain_rows,
907 size_symm);
908 }
909#endif // NDEBUG
910
911 const auto v = getVolume();
912 auto t_w = getFTensor0IntegrationWeight();
913 const int nb_base_functions = data.getN().size2();
914 auto t_row_base_fun = data.getFTensor0N();
915
916 auto t_obj_dstrain =
917 getFTensor2SymmetricFromMat<SPACE_DIM>(topoData->objDStrainAtPts);
918
919 FTENSOR_INDEX(SPACE_DIM, i);
920 FTENSOR_INDEX(SPACE_DIM, j);
921 FTensor::Index<'L', size_symm> L;
922 auto t_L = symm_L_tensor();
923
924 auto get_ftensor1 = [](auto &v) {
925 return FTensor::Tensor1<FTensor::PackPtr<double *, size_symm>, size_symm>(
926 &v[0], &v[1], &v[2], &v[3], &v[4], &v[5]);
927 };
928
929 for (int gg = 0; gg != nb_integration_pts; ++gg) {
930 const double a = v * t_w;
931 auto t_nf = get_ftensor1(nF);
932
933 FTensor::Tensor1<double, size_symm> t_obj_dU;
934 t_obj_dU(L) = t_L(i, j, L) * t_obj_dstrain(i, j);
935
936 int bb = 0;
937 for (; bb != nb_dofs / size_symm; ++bb) {
938 t_nf(L) += a * t_row_base_fun * t_obj_dU(L);
939 ++t_nf;
940 ++t_row_base_fun;
941 }
942 for (; bb != nb_base_functions; ++bb)
943 ++t_row_base_fun;
944
945 ++t_w;
946 ++t_obj_dstrain;
947 }
948
949 MoFEMFunctionReturn(0);
950}
951
952struct OpSensitivityInteriorGradient
953 : public OpAssembleVolumeTopologicalDerivativeImpl {
954 OpSensitivityInteriorGradient(
955 const std::string field_name,
956 boost::shared_ptr<DataAtIntegrationPts> data_ptr,
957 boost::shared_ptr<TopologicalData> topo_ptr,
958 SmartPetscObj<Vec> assemble_vec, const double alpha,
959 const double rho, const double alpha_omega = 0)
960 : OpAssembleVolumeTopologicalDerivativeImpl(field_name, data_ptr,
961 topo_ptr, nullptr,
962 assemble_vec, Tag()),
963 alphaW(alpha), alphaRho(rho), alphaOmega(alpha_omega) {}
964
965 MoFEMErrorCode integrate(EntData &data) {
966 MoFEMFunctionBegin;
967
968#ifndef NDEBUG
969 if (!dataAtPts)
970 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
971 "DataAtIntegrationPts pointer is null");
972#endif // NDEBUG
973
974 const int nb_dofs = data.getIndices().size();
975 if (!nb_dofs)
976 MoFEMFunctionReturnHot(0);
977
978 const int nb_integration_pts = data.getN().size1();
979
980#ifndef NDEBUG
981 if (dataAtPts->divPAtPts.size2() != nb_integration_pts)
982 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
983 "divPAtPts columns (%d) != nb integration points (%d)",
984 static_cast<int>(dataAtPts->divPAtPts.size2()),
985 nb_integration_pts);
986 if (dataAtPts->wL2DotAtPts.size2() != nb_integration_pts)
987 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
988 "wL2DotAtPts columns (%d) != nb integration points (%d)",
989 static_cast<int>(dataAtPts->wL2DotAtPts.size2()),
990 nb_integration_pts);
991 if (dataAtPts->varWL2.size2() != nb_integration_pts)
992 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
993 "varWL2 columns (%d) != nb integration points (%d)",
994 static_cast<int>(dataAtPts->varWL2.size2()), nb_integration_pts);
995 if (dataAtPts->varRotAxis.size2() != nb_integration_pts)
996 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
997 "varRotAxis columns (%d) != nb integration points (%d)",
998 static_cast<int>(dataAtPts->varRotAxis.size2()),
999 nb_integration_pts);
1000 if (dataAtPts->varGradRotAxis.size2() != nb_integration_pts)
1001 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1002 "varGradRotAxis columns (%d) != nb integration points (%d)",
1003 static_cast<int>(dataAtPts->varGradRotAxis.size2()),
1004 nb_integration_pts);
1005 if (dataAtPts->varPiola.size2() != nb_integration_pts)
1006 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1007 "varPiola columns (%d) != nb integration points (%d)",
1008 static_cast<int>(dataAtPts->varPiola.size2()),
1009 nb_integration_pts);
1010 if (dataAtPts->varDivPiola.size2() != nb_integration_pts)
1011 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1012 "varDivPiola columns (%d) != nb integration points (%d)",
1013 static_cast<int>(dataAtPts->varDivPiola.size2()),
1014 nb_integration_pts);
1015#endif // NDEBUG
1016
1017 const auto v = getVolume();
1018 auto t_w = getFTensor0IntegrationWeight();
1019 auto t_div_P = getFTensor1FromMat<SPACE_DIM>(dataAtPts->divPAtPts);
1020 auto t_w_l2 = getFTensor1FromMat<SPACE_DIM>(dataAtPts->wL2AtPts);
1021 auto t_s_dot_w = getFTensor1FromMat<SPACE_DIM>(dataAtPts->wL2DotAtPts);
1022 auto t_s_dot_dot_w =
1023 getFTensor1FromMat<SPACE_DIM>(dataAtPts->wL2DotDotAtPts);
1024 auto t_h = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(dataAtPts->hAtPts);
1025 auto t_levi_kirchhoff =
1026 getFTensor1FromMat<SPACE_DIM>(dataAtPts->leviKirchhoffAtPts);
1027 auto t_omega_grad_dot =
1028 getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(dataAtPts->rotAxisGradDotAtPts);
1029 auto t_R = getFTensor2FromMat<3, 3>(dataAtPts->rotMatAtPts);
1030 auto t_u = getFTensor2SymmetricFromMat<3>(dataAtPts->stretchTensorAtPts);
1031
1032 auto t_var_w = getFTensor1FromMat<SPACE_DIM>(dataAtPts->varWL2);
1033 auto t_var_omega = getFTensor1FromMat<SPACE_DIM>(dataAtPts->varRotAxis);
1034 auto t_var_grad_omega =
1035 getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(dataAtPts->varGradRotAxis);
1036 auto t_var_P = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(dataAtPts->varPiola);
1037 auto t_var_div_P = getFTensor1FromMat<SPACE_DIM>(dataAtPts->varDivPiola);
1038
1039 auto t_det = getFTensor0FromVec(topoData->detJacobianAtPts);
1040 auto t_inv_jac =
1041 getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(topoData->invJacobianAtPtr);
1042
1043 if (dataAtPts->wL2DotDotAtPts.size1() == 0 &&
1044 dataAtPts->wL2DotDotAtPts.size2() != nb_integration_pts) {
1045 dataAtPts->wL2DotDotAtPts.resize(SPACE_DIM, nb_integration_pts);
1046 dataAtPts->wL2DotDotAtPts.clear();
1047 }
1048
1049 const auto piola_scale = dataAtPts->piolaScale;
1050 const auto alpha_w = alphaW / piola_scale;
1051 const auto alpha_rho = alphaRho / piola_scale;
1052
1053 const int nb_base_functions = data.getN().size2();
1054 auto t_base_diff = data.getFTensor1DiffN<3>();
1055
1056 auto get_ftensor1 = [](auto &v) {
1057 return FTensor::Tensor1<FTensor::PackPtr<double *, SPACE_DIM>, SPACE_DIM>(
1058 &v[0], &v[1], &v[2]);
1059 };
1060
1061 auto next = [&]() {
1062 ++t_w;
1063 ++t_div_P;
1064 ++t_w_l2;
1065 ++t_s_dot_w;
1066 ++t_s_dot_dot_w;
1067 ++t_h;
1068 ++t_levi_kirchhoff;
1069 ++t_omega_grad_dot;
1070 ++t_R;
1071 ++t_u;
1072
1073 ++t_var_w;
1074 ++t_var_omega;
1075 ++t_var_grad_omega;
1076 ++t_var_P;
1077 ++t_var_div_P;
1078
1079 ++t_det;
1080 ++t_inv_jac;
1081 };
1082
1083 FTENSOR_INDEX(SPACE_DIM, i);
1084 FTENSOR_INDEX(SPACE_DIM, j);
1085 FTENSOR_INDEX(SPACE_DIM, k);
1086 FTENSOR_INDEX(SPACE_DIM, l);
1087 FTENSOR_INDEX(SPACE_DIM, m);
1088 constexpr auto t_kd = FTensor::Kronecker_Delta<double>();
1089
1090 for (int gg = 0; gg != nb_integration_pts; ++gg) {
1091
1092 // Calculate the variation of the gradient due to geometry change
1093 FTensor::Tensor2<double, SPACE_DIM, SPACE_DIM> t_cof;
1094 t_cof(i, j) = t_det * t_inv_jac(j, i);
1095
1096 auto t_nf = get_ftensor1(nF);
1097 int bb = 0;
1098 for (; bb != nb_dofs / SPACE_DIM; ++bb) {
1099
1100 FTensor::Tensor1<double, SPACE_DIM> t_div_base;
1101 t_div_base(i) = -(1 / t_det) * (t_inv_jac(j, i) * t_base_diff(j));
1102
1103 // OpSpatialEquilibrium
1104 t_nf(i) += (t_w * v) *
1105 (t_var_w(k) * (-t_div_P(k) + alpha_w * t_s_dot_w(k) +
1106 alpha_rho * t_s_dot_dot_w(k))) *
1107 t_cof(i, j) * t_base_diff(j);
1108 t_nf(i) += (t_w * v) * (-(t_var_w(k) * t_div_P(k))) * t_div_base(i);
1109
1110 // OpSpatialRotation
1111 t_nf(i) += (t_w * v) * (t_var_omega(k) * (-t_levi_kirchhoff(k))) *
1112 t_cof(i, j) * t_base_diff(j);
1113 #ifndef NDEBUG
1114 // need to add implementaion of omege terms
1115 if (alphaOmega) {
1116 SETERRQ(
1117 PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
1118 "OpSensitivity_dX with alpha_omega != 0 is not implemented yet");
1119 }
1120 #endif
1121
1122 // OpSpatialConsistencyP
1123 if (EshelbianCore::gradApproximator == NO_H1_CONFIGURATION) {
1124 t_nf(i) -= (t_w * v) * (t_var_P(i, k) * (t_R(i, l) * t_u(l, k)) / 2) *
1125 t_cof(i, j) * t_base_diff(j);
1126 t_nf(i) -= (t_w * v) * (t_var_P(i, l) * (t_R(i, k) * t_u(l, k)) / 2) *
1127 t_cof(i, j) * t_base_diff(j);
1128 t_nf(i) += (t_w * v) * (t_var_P(i, j) * t_kd(i, j)) * t_cof(i, j) *
1129 t_base_diff(j);
1130 } else {
1131 FTensor::Tensor2<double, SPACE_DIM, SPACE_DIM> t_residuum_P;
1132 t_residuum_P(k, m) = t_h(k, m) - t_kd(k, m);
1133 t_nf(i) += (t_w * v) * (t_var_P(k, m) * (-t_residuum_P(k, m))) *
1134 t_cof(i, j) * t_base_diff(j);
1135 }
1136
1137 // OpSpatialConsistencyDivTerm
1138 t_nf(i) += (t_w * v) * (t_var_div_P(k) * (-t_w_l2(k))) * t_cof(i, j) *
1139 t_base_diff(j);
1140 t_nf(i) += (t_w * v) * (t_var_div_P(k) * (-t_w_l2(k))) * t_div_base(i);
1141
1142 ++t_nf;
1143 ++t_base_diff;
1144 }
1145 for (; bb != nb_base_functions; ++bb)
1146 ++t_base_diff;
1147
1148 next();
1149 }
1150
1151 MoFEMFunctionReturn(0);
1152 }
1153
1154private:
1155 const double alphaW;
1156 const double alphaRho;
1157 const double alphaOmega;
1158};
1159
1160} // namespace EshelbianPlasticity
FTENSOR_INDEX
#define FTENSOR_INDEX(DIM, I)
Definition Templates.hpp:2110
a
constexpr double a
Definition approx_sphere.cpp:30
SPACE_DIM
constexpr int SPACE_DIM
Definition child_and_parent.cpp:17
FTensor::Kronecker_Delta
Kronecker Delta class.
Definition Kronecker_Delta.hpp:15
FTensor::Tensor2_symmetric
Definition Tensor2_symmetric_value.hpp:14
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition definitions.h:463
NOSPACE
@ NOSPACE
Definition definitions.h:83
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition definitions.h:359
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition definitions.h:31
MOFEM_NOT_IMPLEMENTED
@ MOFEM_NOT_IMPLEMENTED
Definition definitions.h:32
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition definitions.h:432
CHKERR
#define CHKERR
Inline error check.
Definition definitions.h:551
t_kd
constexpr auto t_kd
Definition free_surface.cpp:163
i
FTensor::Index< 'i', SPACE_DIM > i
Definition hcurl_divergence_operator_2d.cpp:27
v
const double v
phase velocity of light in medium (cm/ns)
Definition initial_diffusion.cpp:40
l
FTensor::Index< 'l', 3 > l
Definition matrix_function.cpp:21
j
FTensor::Index< 'j', 3 > j
Definition matrix_function.cpp:19
k
FTensor::Index< 'k', 3 > k
Definition matrix_function.cpp:20
EshelbianPlasticity::size_symm
static constexpr auto size_symm
Definition EshelbianAux.hpp:41
A
constexpr AssemblyType A
Definition operators_tests.cpp:30
field_name
constexpr auto field_name
Definition poisson_2d_homogeneous.cpp:13
m
FTensor::Index< 'm', 3 > m
Definition shallow_wave.cpp:80
EshelbianCore::gradApproximator
static enum RotSelector gradApproximator
Definition EshelbianCore.hpp:28
EshelbianPlasticity::OpAssembleBrokenFaceTopologicalDerivativeImplBase::OpAssembleBrokenFaceTopologicalDerivativeImplBase
OpAssembleBrokenFaceTopologicalDerivativeImplBase(boost::shared_ptr< std::vector< BrokenBaseSideData > > broken_base_side_data, boost::shared_ptr< TopologicalData > topo_ptr, boost::shared_ptr< double > J_ptr, SmartPetscObj< Vec > assemble_vec, Tag topo_tag, boost::shared_ptr< Range > ents_ptr=nullptr)
Definition EshelbianTopologicalDerivativeOperators.cpp:103
EshelbianPlasticity::OpAssembleTopologicalObjectiveDerivativeImplBase::OpAssembleTopologicalObjectiveDerivativeImplBase
OpAssembleTopologicalObjectiveDerivativeImplBase(const std::string &field_name, boost::shared_ptr< DataAtIntegrationPts > data_ptr, boost::shared_ptr< TopologicalData > topo_ptr, boost::shared_ptr< double > J_ptr, SmartPetscObj< Vec > assemble_vec, Tag topo_tag)
Definition EshelbianTopologicalDerivativeOperators.cpp:36
EshelbianPlasticity::OpJ_dBubbleImpl
Definition EshelbianTopologicalDerivativeOperators.cpp:156
EshelbianPlasticity::OpJ_dBubbleImpl::integrate
MoFEMErrorCode integrate(EntData &data)
Definition EshelbianTopologicalDerivativeOperators.cpp:461
EshelbianPlasticity::OpJ_dBubble_no_streachImpl
Definition EshelbianTopologicalDerivativeOperators.cpp:194
EshelbianPlasticity::OpJ_dBubble_no_streachImpl::integrate
MoFEMErrorCode integrate(EntData &data)
Definition EshelbianTopologicalDerivativeOperators.cpp:797
EshelbianPlasticity::OpJ_dPImpl
Definition EshelbianTopologicalDerivativeOperators.cpp:149
EshelbianPlasticity::OpJ_dPImpl::integrate
MoFEMErrorCode integrate(EntData &data)
Definition EshelbianTopologicalDerivativeOperators.cpp:395
EshelbianPlasticity::OpJ_dP_no_streachImpl
Definition EshelbianTopologicalDerivativeOperators.cpp:186
EshelbianPlasticity::OpJ_dP_no_streachImpl::integrate
MoFEMErrorCode integrate(EntData &data)
Definition EshelbianTopologicalDerivativeOperators.cpp:712
EshelbianPlasticity::OpJ_dUImpl
Definition EshelbianTopologicalDerivativeOperators.cpp:201
EshelbianPlasticity::OpJ_dUImpl::integrate
MoFEMErrorCode integrate(EntData &data)
Definition EshelbianTopologicalDerivativeOperators.cpp:881
EshelbianPlasticity::OpJ_dwImpl
Definition EshelbianTopologicalDerivativeOperators.cpp:163
EshelbianPlasticity::OpJ_dwImpl::integrate
MoFEMErrorCode integrate(EntData &data)
Definition EshelbianTopologicalDerivativeOperators.cpp:525
EshelbianPlasticity::OpSensitivityInteriorGradient::OpSensitivityInteriorGradient
OpSensitivityInteriorGradient(const std::string field_name, boost::shared_ptr< DataAtIntegrationPts > data_ptr, boost::shared_ptr< TopologicalData > topo_ptr, SmartPetscObj< Vec > assemble_vec, const double alpha, const double rho, const double alpha_omega=0)
Definition EshelbianTopologicalDerivativeOperators.cpp:954
EshelbianPlasticity::OpTopologicalObjectivePythonImpl::OpTopologicalObjectivePythonImpl
OpTopologicalObjectivePythonImpl(boost::shared_ptr< DataAtIntegrationPts > data_at_pts_ptr, boost::shared_ptr< TopologicalData > topo_p, boost::shared_ptr< ObjectiveFunctionData > python_ptr)
Definition EshelbianTopologicalDerivativeOperators.cpp:219
EshelbianPlasticity::OpTopologicalObjectivePythonImpl::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Definition EshelbianTopologicalDerivativeOperators.cpp:234
EshelbianPlasticity::dJ_dTractionImpl
Definition EshelbianTopologicalDerivativeOperators.cpp:178
EshelbianPlasticity::dJ_dTractionImpl::iNtegrate
MoFEMErrorCode iNtegrate(EntData &data)
Definition EshelbianTopologicalDerivativeOperators.cpp:652
EshelbianPlasticity::dJ_duGammaImpl
Definition EshelbianTopologicalDerivativeOperators.cpp:170
EshelbianPlasticity::dJ_duGammaImpl::integrate
MoFEMErrorCode integrate(EntData &data)
Definition EshelbianTopologicalDerivativeOperators.cpp:589
MoFEM::EntitiesFieldData::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition EntitiesFieldData.hpp:152
MoFEM::EntitiesFieldData::EntData::getFTensor1DiffN
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
Definition EntitiesFieldData.cpp:533
MoFEM::EntitiesFieldData::EntData::getFTensor2N
FTensor::Tensor2< FTensor::PackPtr< double *, Tensor_Dim0 *Tensor_Dim1 >, Tensor_Dim0, Tensor_Dim1 > getFTensor2N(FieldApproximationBase base)
Get base functions for Hdiv/Hcurl spaces.
Definition EntitiesFieldData.cpp:769
MoFEM::EntitiesFieldData::EntData::getFTensor0N
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
Definition EntitiesFieldData.hpp:1692
MoFEM::EntitiesFieldData::EntData::getFieldEntities
const VectorFieldEntities & getFieldEntities() const
Get field entities (const version)
Definition EntitiesFieldData.hpp:1450
MoFEM::EntitiesFieldData::EntData::getN
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
Definition EntitiesFieldData.hpp:1486
MoFEM::EntitiesFieldData::EntData::getFTensor1N
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1N(FieldApproximationBase base)
Get base functions for Hdiv/Hcurl spaces.
Definition EntitiesFieldData.cpp:647
MoFEM::EntitiesFieldData::EntData::getIndices
const VectorInt & getIndices() const
Get global indices of degrees of freedom on entity.
Definition EntitiesFieldData.hpp:1382
MoFEM::ForcesAndSourcesCore::UserDataOperator::getFTensor0IntegrationWeight
auto getFTensor0IntegrationWeight()
Get integration weights.
Definition ForcesAndSourcesCore.hpp:1354
OpAssembleBasic< VolUserDataOperator >::nF
VectorDouble nF
local right hand side vector
Definition EshelbianOperators.hpp:62
OpAssembleBasic< VolUserDataOperator >::dataAtPts
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
data at integration pts
Definition EshelbianOperators.hpp:44
rho
double rho
Definition plastic.cpp:144
Generated by Doxygen 1.9.8 and hosted at University of Glasgow