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