v0.15.0
Loading...
Searching...
No Matches
adjoint.cpp
Go to the documentation of this file.
1/**
2 * @file adjoint.cpp
3 * @brief Topology optimization using adjoint method with Python objective functions
4 * @details This tutorial demonstrates:
5 * - Adjoint-based sensitivity analysis for structural topology optimization
6 * - Integration with Python for flexible objective function definition
7 * - Higher-order finite element analysis with geometry fields
8 * - Automatic differentiation using MoFEM's adjoint capabilities
9 * @author Anonymous author(s) committing under MIT license
10 * @example mofem/tutorials/vec-7/adjoint.cpp
11 *
12 */
13
14#include <boost/python.hpp>
15#include <boost/python/def.hpp>
16#include <boost/python/numpy.hpp>
17namespace bp = boost::python;
18namespace np = boost::python::numpy;
19
20#include <MoFEM.hpp>
21
22using namespace MoFEM;
23
24//! [Constants and material properties]
25constexpr int BASE_DIM = 1; ///< Dimension of the base functions
26
27//! [Define dimension]
28constexpr int SPACE_DIM =
29 EXECUTABLE_DIMENSION; ///< Space dimension of problem (2D or 3D), set at compile time
30
31//! [Define dimension]
32constexpr AssemblyType A = AssemblyType::PETSC; ///< Use PETSc for matrix/vector assembly
33constexpr IntegrationType I =
34 IntegrationType::GAUSS; ///< Use Gauss quadrature for integration
35
36//! [Material properties for linear elasticity]
37constexpr double young_modulus = 1; ///< Young's modulus E
38constexpr double poisson_ratio = 0.3; ///< Poisson's ratio ν
39constexpr double bulk_modulus_K = young_modulus / (3 * (1 - 2 * poisson_ratio)); ///< Bulk modulus K = E/(3(1-2ν))
40constexpr double shear_modulus_G = young_modulus / (2 * (1 + poisson_ratio)); ///< Shear modulus G = E/(2(1+ν))
41
42PetscBool is_plane_strain = PETSC_FALSE; ///< Flag for plane strain vs plane stress in 2D
43//! [Constants and material properties]
44
45//! [Define finite element types and operators]
46using EntData = EntitiesFieldData::EntData; ///< Entity data for field operations
49using DomainEleOp = DomainEle::UserDataOperator; ///< Domain element operators
50using BoundaryEleOp = BoundaryEle::UserDataOperator; ///< Boundary element operators
51//! [Define finite element types and operators]
52
53//! [Boundary condition types]
54struct DomainBCs {}; ///< Domain boundary conditions marker
55struct BoundaryBCs {}; ///< Boundary conditions marker
56
57//! [Natural boundary condition operators]
64
65template <int DIM> struct PostProcEleByDim;
66
67template <> struct PostProcEleByDim<2> {
71};
72
73template <> struct PostProcEleByDim<3> {
77};
78
82
83//! [Adjoint operator declarations]
84/// Forward declaration of adjoint test operator for geometric sensitivity
85template <int SPACE_DIM, IntegrationType I, typename OpBase>
87
88/// Forward declaration of operator for gradient times symmetric tensor operations
89template <int SPACE_DIM, IntegrationType I, typename OpBase>
91//! [Adjoint operator declarations]
92
93constexpr double young_modulus = 1;
94constexpr double poisson_ratio = 0.3;
95constexpr double bulk_modulus_K = young_modulus / (3 * (1 - 2 * poisson_ratio));
96constexpr double shear_modulus_G = young_modulus / (2 * (1 + poisson_ratio));
97
98PetscBool is_plane_strain = PETSC_FALSE;
99
100#include <ElasticSpring.hpp>
101#include <FluidLevel.hpp>
102#include <CalculateTraction.hpp>
103#include <NaturalDomainBC.hpp>
104#include <NaturalBoundaryBC.hpp>
105#include <HookeOps.hpp>
106
107/**
108 * @brief Abstract interface for Python-defined objective functions
109 *
110 * This interface allows users to define custom objective functions in Python
111 * that can be called from C++ during the optimization process. The objective
112 * function and its gradients are evaluated at integration points using:
113 * - Coordinates of integration points
114 * - Displacement field values
115 * - Stress and strain tensors
116 *
117 * Key methods:
118 * - evalObjectiveFunction: Compute f(x,u,σ,ε)
119 * - evalObjectiveGradientStress: Compute ∂f/∂σ
120 * - evalObjectiveGradientStrain: Compute ∂f/∂ε
121 * - numberOfModes: Return number of topology modes for optimization
122 * - blockModes: Define spatial topology modes for shape optimization
123 */
125
126 /// Evaluate objective function f(coords, u, stress, strain)
127 virtual MoFEMErrorCode
129 boost::shared_ptr<MatrixDouble> u_ptr,
130 boost::shared_ptr<MatrixDouble> stress_ptr,
131 boost::shared_ptr<MatrixDouble> strain_ptr,
132 boost::shared_ptr<VectorDouble> o_ptr) = 0;
133
134 /// Evaluate gradient of objective function w.r.t. stress: ∂f/∂σ
135 virtual MoFEMErrorCode
137 boost::shared_ptr<MatrixDouble> u_ptr,
138 boost::shared_ptr<MatrixDouble> stress_ptr,
139 boost::shared_ptr<MatrixDouble> strain_ptr,
140 boost::shared_ptr<MatrixDouble> o_ptr) = 0;
141
142 /// Evaluate gradient of objective function w.r.t. strain: ∂f/∂ε
143 virtual MoFEMErrorCode
145 boost::shared_ptr<MatrixDouble> u_ptr,
146 boost::shared_ptr<MatrixDouble> stress_ptr,
147 boost::shared_ptr<MatrixDouble> strain_ptr,
148 boost::shared_ptr<MatrixDouble> o_ptr) = 0;
149
150 /// Return number of topology optimization modes for given block
151 virtual MoFEMErrorCode numberOfModes(int block_id, int &modes) = 0;
152
153 /// Define spatial modes for topology optimization
154 virtual MoFEMErrorCode blockModes(int block_id, MatrixDouble &coords,
155 std::array<double, 3> &centroid,
156 std::array<double, 6> &bbox,
157 MatrixDouble &o_ptr) = 0;
158
159 virtual ~ObjectiveFunctionData() = default;
160};
161
162boost::shared_ptr<ObjectiveFunctionData>
163create_python_objective_function(std::string py_file);
164
166 const std::string block_name, int dim);
167
168MoFEMErrorCode save_range(moab::Interface &moab, const std::string name,
169 const Range r);
170/**
171 * @brief Main class for topology optimization using adjoint method
172 *
173 * This class implements a complete topology optimization workflow:
174 * 1. Setup finite element problems for elasticity and adjoint analysis
175 * 2. Define topology modes for design parametrization
176 * 3. Use TAO optimization library to minimize Python-defined objective function
177 * 4. Compute sensitivities using adjoint method for efficient gradient calculation
178 *
179 * The optimization process alternates between:
180 * - Forward analysis: solve elastic problem for current geometry
181 * - Adjoint analysis: solve adjoint problem for objective function gradients
182 * - Geometry update: modify design using optimization algorithm (L-BFGS)
183 */
184struct Example {
185
186 Example(MoFEM::Interface &m_field) : mField(m_field) {}
187
188 /// Main driver function for the optimization process
190
191private:
192 MoFEM::Interface &mField; ///< Reference to MoFEM interface
193
194 boost::shared_ptr<MatrixDouble> vectorFieldPtr = nullptr; ///< Field values at evaluation points
195
196 // Problem setup methods
197 MoFEMErrorCode readMesh(); ///< Read mesh from file and setup meshsets
198 MoFEMErrorCode setupProblem(); ///< Setup fields, approximation spaces and DOFs
199 MoFEMErrorCode setupAdJoint(); ///< Setup adjoint fields and finite elements
200 MoFEMErrorCode boundaryCondition(); ///< Apply essential boundary conditions
201 MoFEMErrorCode topologyModes(); ///< Compute topology optimization modes
202 MoFEMErrorCode assembleSystem(); ///< Setup operators in finite element pipeline
203
204 // Analysis methods
205 MoFEMErrorCode solveElastic(); ///< Solve forward elastic problem
206 MoFEMErrorCode postprocessElastic(int iter, SmartPetscObj<Vec> adjoint_vector = nullptr); ///< Post-process and output results
207
208 /// Calculate objective function gradient using adjoint method
209 MoFEMErrorCode calculateGradient(PetscReal *objective_function_value,
210 Vec objective_function_gradient,
211 Vec adjoint_vector);
212
213 // Problem configuration
214 FieldApproximationBase base; ///< Choice of finite element basis functions
215 int fieldOrder = 2; ///< Polynomial order for approximation
216
217 // Solver and data management
218 SmartPetscObj<KSP> kspElastic; ///< Linear solver for elastic problem
219 SmartPetscObj<DM> adjointDM; ///< Data manager for adjoint problem
220 boost::shared_ptr<ObjectiveFunctionData> pythonPtr; ///< Interface to Python objective function
221
222 // Topology optimization data
223 std::vector<SmartPetscObj<Vec>> modeVecs; ///< Topology mode vectors (design variables)
224 std::vector<std::array<double, 3>> modeCentroids; ///< Centroids of optimization blocks
225 std::vector<std::array<double, 6>> modeBBoxes; ///< Bounding boxes of optimization blocks
226 SmartPetscObj<Vec> initialGeometry; ///< Initial geometry field
227};
228
229//! [Run topology optimization problem]
230/**
231 * @brief Main driver for topology optimization using adjoint sensitivity analysis
232 *
233 * This function orchestrates the complete topology optimization workflow:
234 * 1. Initialize Python objective function interface
235 * 2. Setup finite element problems (forward and adjoint)
236 * 3. Compute initial elastic solution
237 * 4. Generate topology optimization modes
238 * 5. Run TAO optimization loop with adjoint-based gradients
239 *
240 * The optimization uses TAO (Toolkit for Advanced Optimization) with L-BFGS
241 * algorithm. At each iteration:
242 * - Update geometry based on current design variables
243 * - Solve forward elastic problem
244 * - Compute objective function and gradients using adjoint method
245 * - Post-process results
246 *
247 * @return MoFEMErrorCode Success or error code
248 */
251
252 // Read objective function from Python file
253 char objective_function_file_name[255] = "objective_function.py";
255 PETSC_NULLPTR, PETSC_NULLPTR, "-objective_function",
256 objective_function_file_name, 255, PETSC_NULLPTR);
257
258 // Verify that the Python objective function file exists
259 auto file_exists = [](std::string myfile) {
260 std::ifstream file(myfile.c_str());
261 if (file) {
262 return true;
263 }
264 return false;
265 };
266 if (!file_exists(objective_function_file_name)) {
267 MOFEM_LOG("WORLD", Sev::error) << "Objective function file NOT found: "
268 << objective_function_file_name;
269 CHK_THROW_MESSAGE(MOFEM_NOT_FOUND, "file NOT found");
270 }
271
272 // Create Python interface for objective function
273 pythonPtr = create_python_objective_function(objective_function_file_name);
274
275 // Setup finite element problems
276 CHKERR readMesh(); // Read mesh and meshsets
277 CHKERR setupProblem(); // Setup displacement field and geometry field
278 CHKERR setupAdJoint(); // Setup adjoint field for sensitivity analysis
279 CHKERR boundaryCondition(); // Apply essential boundary conditions
280 CHKERR assembleSystem(); // Setup finite element operators
281
282 // Create linear solver for elastic problem
283 auto pip = mField.getInterface<PipelineManager>();
284 kspElastic = pip->createKSP();
285 CHKERR KSPSetFromOptions(kspElastic);
286
287 // Solve initial elastic problem
289 CHKERR postprocessElastic(-1); // Post-process initial solution
290
291 auto create_vec_modes = [&](auto block_name) {
293 auto mesh_mng = mField.getInterface<MeshsetsManager>();
294 auto bcs = mesh_mng->getCubitMeshsetPtr(
295
296 std::regex((boost::format("%s(.*)") % block_name).str())
297
298 );
299
300 int nb_total_modes = 0;
301 for (auto &bc : bcs) {
302 auto id = bc->getMeshsetId();
303 int nb_modes;
304 CHKERR pythonPtr->numberOfModes(id, nb_modes);
305 nb_total_modes += nb_modes;
306 }
307
308 MOFEM_LOG("WORLD", Sev::inform)
309 << "Total number of modes to apply: " << nb_total_modes;
310
311 modeVecs.resize(nb_total_modes);
312
314 };
315
316 auto get_modes_bounding_boxes = [&](auto block_name) {
318 auto mesh_mng = mField.getInterface<MeshsetsManager>();
319 auto bcs = mesh_mng->getCubitMeshsetPtr(
320
321 std::regex((boost::format("%s(.*)") % block_name).str())
322
323 );
324
325 for (auto &bc : bcs) {
326 auto meshset = bc->getMeshset();
327 Range ents;
328 CHKERR mField.get_moab().get_entities_by_handle(meshset, ents, true);
329 Range verts;
330 CHKERR mField.get_moab().get_connectivity(ents, verts, false);
331 std::vector<double> x(verts.size());
332 std::vector<double> y(verts.size());
333 std::vector<double> z(verts.size());
334 CHKERR mField.get_moab().get_coords(verts, x.data(), y.data(), z.data());
335 std::array<double, 3> centroid = {0, 0, 0};
336 for (int i = 0; i != verts.size(); ++i) {
337 centroid[0] += x[i];
338 centroid[1] += y[i];
339 centroid[2] += z[i];
340 }
341 MPI_Allreduce(MPI_IN_PLACE, centroid.data(), 3, MPI_DOUBLE, MPI_SUM,
342 mField.get_comm());
343 int nb_vertex = verts.size();
344 MPI_Allreduce(MPI_IN_PLACE, &nb_vertex, 1, MPI_INT, MPI_SUM,
345 mField.get_comm());
346 if (nb_vertex) {
347 centroid[0] /= nb_vertex;
348 centroid[1] /= nb_vertex;
349 centroid[2] /= nb_vertex;
350 }
351 std::array<double, 6> bbox = {centroid[0], centroid[1], centroid[2],
352 centroid[0], centroid[1], centroid[2]};
353 for (int i = 0; i != verts.size(); ++i) {
354 bbox[0] = std::min(bbox[0], x[i]);
355 bbox[1] = std::min(bbox[1], y[i]);
356 bbox[2] = std::min(bbox[2], z[i]);
357 bbox[3] = std::max(bbox[3], x[i]);
358 bbox[4] = std::max(bbox[4], y[i]);
359 bbox[5] = std::max(bbox[5], z[i]);
360 }
361 MPI_Allreduce(MPI_IN_PLACE, &bbox[0], 3, MPI_DOUBLE, MPI_MIN,
362 mField.get_comm());
363 MPI_Allreduce(MPI_IN_PLACE, &bbox[3], 3, MPI_DOUBLE, MPI_MAX,
364 mField.get_comm());
365
366 MOFEM_LOG("WORLD", Sev::inform)
367 << "Block: " << bc->getName() << " centroid: " << centroid[0] << " "
368 << centroid[1] << " " << centroid[2];
369 MOFEM_LOG("WORLD", Sev::inform)
370 << "Block: " << bc->getName() << " bbox: " << bbox[0] << " "
371 << bbox[1] << " " << bbox[2] << " " << bbox[3] << " " << bbox[4]
372 << " " << bbox[5];
373
374 modeCentroids.push_back(centroid);
375 modeBBoxes.push_back(bbox);
376 }
377
379 };
380
381 auto eval_objective_and_gradient = [](Tao tao, Vec x, PetscReal *f, Vec g,
382 void *ctx) -> PetscErrorCode {
384 auto *ex_ptr = (Example *)ctx;
385
386 int iter;
387 CHKERR TaoGetIterationNumber(tao, &iter);
388
389 auto set_geometry = [&](Vec x) {
391
392 VecScatter ctx;
393 Vec x_local;
394 CHKERR VecScatterCreateToAll(x, &ctx, &x_local);
395 // scatter as many times as you need
396 CHKERR VecScatterBegin(ctx, x, x_local, INSERT_VALUES, SCATTER_FORWARD);
397 CHKERR VecScatterEnd(ctx, x, x_local, INSERT_VALUES, SCATTER_FORWARD);
398 // destroy scatter context and local vector when no longer needed
399 CHKERR VecScatterDestroy(&ctx);
400
401 auto current_geometry = vectorDuplicate(ex_ptr->initialGeometry);
402 CHKERR VecCopy(ex_ptr->initialGeometry, current_geometry);
403 const double *a;
404 CHKERR VecGetArrayRead(x_local, &a);
405 const double *coeff = a;
406 for (auto &mode_vec : ex_ptr->modeVecs) {
407 MOFEM_LOG("WORLD", Sev::verbose)
408 << "Adding mode with coeff: " << *coeff;
409 CHKERR VecAXPY(current_geometry, (*coeff), mode_vec);
410 ++coeff;
411 }
412 CHKERR VecRestoreArrayRead(x_local, &a);
413 CHKERR VecGhostUpdateBegin(current_geometry, INSERT_VALUES,
414 SCATTER_FORWARD);
415 CHKERR VecGhostUpdateEnd(current_geometry, INSERT_VALUES,
416 SCATTER_FORWARD);
417 CHKERR ex_ptr->mField.getInterface<VecManager>()
418 ->setOtherLocalGhostVector("ADJOINT", "ADJOINT_FIELD", "GEOMETRY",
419 RowColData::ROW, current_geometry,
420 INSERT_VALUES, SCATTER_REVERSE);
421
422 CHKERR VecDestroy(&x_local);
424 };
425
426 CHKERR set_geometry(x);
427 CHKERR KSPReset(ex_ptr->kspElastic);
428 CHKERR ex_ptr->solveElastic();
429 auto simple = ex_ptr->mField.getInterface<Simple>();
430 auto adjoint_vector = createDMVector(simple->getDM());
431 CHKERR ex_ptr->calculateGradient(f, g, adjoint_vector);
432 CHKERR ex_ptr->postprocessElastic(iter, adjoint_vector);
433
435 };
436
437 // Create topology optimization modes and setup TAO solver
438 CHKERR create_vec_modes("OPTIMISE");
439 CHKERR get_modes_bounding_boxes("OPTIMISE");
440
441 /**
442 * Setup TAO (Toolkit for Advanced Optimization) solver for topology optimization
443 *
444 * TAO provides various optimization algorithms. Here we use TAOLMVM (Limited Memory
445 * Variable Metric) which is a quasi-Newton method suitable for unconstrained
446 * optimization with gradient information.
447 *
448 * The optimization variables are coefficients for the topology modes, and
449 * gradients are computed using the adjoint method for efficiency.
450 */
451 auto tao = createTao(mField.get_comm());
452 CHKERR TaoSetType(tao, TAOLMVM);
453
454 auto rank = mField.get_comm_rank();
455 auto g = createVectorMPI(mField.get_comm(), (!rank) ? modeVecs.size() : 0,
456 PETSC_DECIDE);
457 CHKERR TaoSetObjectiveAndGradient(tao, g, eval_objective_and_gradient,
458 (void *)this);
459
460 // Store initial geometry for reference during optimization
462 CHKERR mField.getInterface<VecManager>()->setOtherLocalGhostVector(
463 "ADJOINT", "ADJOINT_FIELD", "GEOMETRY", RowColData::ROW, initialGeometry,
464 INSERT_VALUES, SCATTER_FORWARD);
465
466 // Setup topology optimization modes
467 /**
468 * Generate modes for topology optimization design parameterization
469 *
470 * These modes represent perturbations to the geometry that can be used
471 * as design variables in topology optimization. The modes are defined
472 * through the Python interface and provide spatial basis functions
473 * for shape modifications.
474 */
476
477 // Initialize optimization variables and run solver
478 /**
479 * Start optimization with zero initial guess for design variables
480 *
481 * The TAO solver will iteratively:
482 * 1. Evaluate objective function at current design point
483 * 2. Compute gradients using adjoint sensitivity analysis
484 * 3. Update design variables using L-BFGS algorithm
485 * 4. Check convergence criteria
486 *
487 * Each function evaluation involves solving the forward elastic problem
488 * and the adjoint problem to compute sensitivities efficiently.
489 */
490 auto x0 = vectorDuplicate(g);
491
492 CHKERR VecSet(x0, 0.0);
493 CHKERR TaoSetSolution(tao, x0);
494 CHKERR TaoSetFromOptions(tao);
495 CHKERR TaoSolve(tao);
496
497 // Optimization complete - results available in solution vectors
498 MOFEM_LOG("WORLD", Sev::inform) << "Topology optimization completed";
499
501}
502//! [Run problem]
503
504//! [Read mesh]
505/**
506 * @brief Read mesh from file and setup material/boundary condition meshsets
507 *
508 * This function loads the finite element mesh from file and processes
509 * associated meshsets that define material properties and boundary conditions.
510 * The mesh is typically generated using CUBIT and exported in .h5m format.
511 *
512 * Meshsets are used to group elements/faces by:
513 * - Material properties (for different material blocks)
514 * - Boundary conditions (for applying loads and constraints)
515 * - Optimization regions (for topology optimization)
516 *
517 * @return MoFEMErrorCode Success or error code
518 */
522 CHKERR simple->getOptions(); // Read command line options
523 CHKERR simple->loadFile(); // Load mesh from file
524 // Add meshsets if config file provided
525 CHKERR mField.getInterface<MeshsetsManager>()->setMeshsetFromFile();
527}
528//! [Read mesh]
529
530//! [Set up problem]
531/**
532 * @brief Setup finite element fields, approximation spaces and degrees of freedom
533 *
534 * This function configures the finite element problem by:
535 * 1. Setting up the displacement field "U" with vector approximation
536 * 2. Setting up the geometry field "GEOMETRY" for mesh deformation
537 * 3. Defining polynomial approximation order and basis functions
538 * 4. Creating degrees of freedom on mesh entities
539 *
540 * The displacement field uses H1 vector space for standard elasticity.
541 * The geometry field allows mesh modification during topology optimization.
542 * Different basis functions (Ainsworth-Legendre vs Demkowicz) can be selected.
543 *
544 * @return MoFEMErrorCode Success or error code
545 */
549
550 // Select basis functions for finite element approximation
551 enum bases { AINSWORTH, DEMKOWICZ, LASBASETOPT };
552 const char *list_bases[LASBASETOPT] = {"ainsworth", "demkowicz"};
553 PetscInt choice_base_value = AINSWORTH;
554 CHKERR PetscOptionsGetEList(PETSC_NULLPTR, NULL, "-base", list_bases,
555 LASBASETOPT, &choice_base_value, PETSC_NULLPTR);
556
557 switch (choice_base_value) {
558 case AINSWORTH:
560 MOFEM_LOG("WORLD", Sev::inform)
561 << "Set AINSWORTH_LEGENDRE_BASE for displacements";
562 break;
563 case DEMKOWICZ:
565 MOFEM_LOG("WORLD", Sev::inform)
566 << "Set DEMKOWICZ_JACOBI_BASE for displacements";
567 break;
568 default:
569 base = LASTBASE;
570 break;
571 }
572
573 // Add finite element fields
574 /**
575 * Setup displacement field "U" - the primary unknown in elasticity
576 * This field represents displacement vector at each node/DOF
577 */
578 CHKERR simple->addDomainField("U", H1, base, SPACE_DIM);
579 CHKERR simple->addBoundaryField("U", H1, base, SPACE_DIM);
580 CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-order", &fieldOrder,
581 PETSC_NULLPTR);
582
583 /**
584 * Setup geometry field "GEOMETRY" - used for mesh deformation in optimization
585 * This field stores current nodal coordinates and can be modified
586 * during topology optimization to represent design changes
587 */
588 CHKERR simple->addDataField("GEOMETRY", H1, base, SPACE_DIM);
589
590 // Set polynomial approximation order for both fields
591 CHKERR simple->setFieldOrder("U", fieldOrder);
592 CHKERR simple->setFieldOrder("GEOMETRY", fieldOrder);
593 CHKERR simple->setUp();
594
595 // Project higher-order geometry representation onto geometry field
596 /**
597 * For higher-order elements, this projects the exact geometry
598 * onto the geometry field to maintain curved boundaries accurately
599 */
600 auto project_ho_geometry = [&]() {
601 Projection10NodeCoordsOnField ent_method(mField, "GEOMETRY");
602 return mField.loop_dofs("GEOMETRY", ent_method);
603 };
604 CHKERR project_ho_geometry();
605
606 // Check if plane strain assumption should be used
607 CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-plane_strain",
608 &is_plane_strain, PETSC_NULLPTR);
609
611}
612//! [Set up problem]
613
614//! [Setup adjoint]
615/**
616 * @brief Setup adjoint fields and finite elements for sensitivity analysis
617 *
618 * The adjoint method is used to efficiently compute gradients of the objective
619 * function with respect to design variables. This function sets up:
620 *
621 * 1. ADJOINT_FIELD - stores adjoint variables (Lagrange multipliers)
622 * 2. ADJOINT_DM - data manager for adjoint problem
623 * 3. Adjoint finite elements for domain and boundary
624 *
625 * The adjoint equation is: K^T * λ = ∂f/∂u
626 * where λ are adjoint variables, K is stiffness matrix, f is objective
627 *
628 * @return MoFEMErrorCode Success or error code
629 */
633
634 // Create adjoint data manager and field
635 auto create_adjoint_dm = [&]() {
636 auto adjoint_dm = createDM(mField.get_comm(), "DMMOFEM");
637
638 auto add_field = [&]() {
640 CHKERR mField.add_field("ADJOINT_FIELD", H1, base, SPACE_DIM);
642 "ADJOINT_FIELD");
643 for (auto tt = MBEDGE; tt <= moab::CN::TypeDimensionMap[SPACE_DIM].second;
644 ++tt)
645 CHKERR mField.set_field_order(simple->getMeshset(), tt, "ADJOINT_FIELD",
646 fieldOrder);
647 CHKERR mField.set_field_order(simple->getMeshset(), MBVERTEX,
648 "ADJOINT_FIELD", 1);
651 };
652
653 auto add_adjoint_fe_impl = [&]() {
655 CHKERR mField.add_finite_element("ADJOINT_DOMAIN_FE");
657 "ADJOINT_FIELD");
659 "ADJOINT_FIELD");
661 "ADJOINT_FIELD");
663 "GEOMETRY");
665 simple->getMeshset(), SPACE_DIM, "ADJOINT_DOMAIN_FE");
666 CHKERR mField.build_finite_elements("ADJOINT_DOMAIN_FE");
667
668 CHKERR mField.add_finite_element("ADJOINT_BOUNDARY_FE");
670 "ADJOINT_FIELD");
672 "ADJOINT_FIELD");
674 "ADJOINT_FIELD");
676 "GEOMETRY");
677
678 auto block_name = "OPTIMISE";
679 auto mesh_mng = mField.getInterface<MeshsetsManager>();
680 auto bcs = mesh_mng->getCubitMeshsetPtr(
681
682 std::regex((boost::format("%s(.*)") % block_name).str())
683
684 );
685
686 for (auto bc : bcs) {
688 bc->getMeshset(), SPACE_DIM - 1, "ADJOINT_BOUNDARY_FE");
689 }
690
691 CHKERR mField.build_finite_elements("ADJOINT_BOUNDARY_FE");
692
693 CHKERR mField.build_adjacencies(simple->getBitRefLevel(),
694 simple->getBitRefLevelMask());
695
697 };
698
699 auto set_adjoint_dm_imp = [&]() {
701 CHKERR DMMoFEMCreateMoFEM(adjoint_dm, &mField, "ADJOINT",
702 simple->getBitRefLevel(),
703 simple->getBitRefLevelMask());
704 CHKERR DMMoFEMSetDestroyProblem(adjoint_dm, PETSC_TRUE);
705 CHKERR DMSetFromOptions(adjoint_dm);
706 CHKERR DMMoFEMAddElement(adjoint_dm, "ADJOINT_DOMAIN_FE");
707 CHKERR DMMoFEMAddElement(adjoint_dm, "ADJOINT_BOUNDARY_FE");
708 CHKERR DMMoFEMSetSquareProblem(adjoint_dm, PETSC_TRUE);
709 CHKERR DMMoFEMSetIsPartitioned(adjoint_dm, PETSC_TRUE);
710 mField.getInterface<ProblemsManager>()->buildProblemFromFields =
711 PETSC_TRUE;
712 CHKERR DMSetUp(adjoint_dm);
713 mField.getInterface<ProblemsManager>()->buildProblemFromFields =
714 PETSC_FALSE;
716 };
717
718 CHK_THROW_MESSAGE(add_field(), "add adjoint field");
719 CHK_THROW_MESSAGE(add_adjoint_fe_impl(), "add adjoint fe");
720 CHK_THROW_MESSAGE(set_adjoint_dm_imp(), "set adjoint dm");
721
722 return adjoint_dm;
723 };
724
725 adjointDM = create_adjoint_dm();
726
728}
729//
730
731//! [Boundary condition]
735 auto bc_mng = mField.getInterface<BcManager>();
736
737 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_X",
738 "U", 0, 0);
739 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_Y",
740 "U", 1, 1);
741 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_Z",
742 "U", 2, 2);
743 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(),
744 "REMOVE_ALL", "U", 0, 3);
745 CHKERR bc_mng->pushMarkDOFsOnEntities<DisplacementCubitBcData>(
746 simple->getProblemName(), "U");
747
749}
750//! [Boundary condition]
751
752//! [Adjoint modes]
755
756 auto opt_ents = get_range_from_block(mField, "OPTIMISE", SPACE_DIM - 1);
757 auto subset_dm_bdy = createDM(mField.get_comm(), "DMMOFEM");
758 CHKERR DMMoFEMSetSquareProblem(subset_dm_bdy, PETSC_TRUE);
759 CHKERR DMMoFEMCreateSubDM(subset_dm_bdy, adjointDM, "SUBSET_BDY");
760 CHKERR DMMoFEMAddElement(subset_dm_bdy, "ADJOINT_BOUNDARY_FE");
761 CHKERR DMMoFEMAddSubFieldRow(subset_dm_bdy, "ADJOINT_FIELD",
762 boost::make_shared<Range>(opt_ents));
763 CHKERR DMMoFEMAddSubFieldCol(subset_dm_bdy, "ADJOINT_FIELD",
764 boost::make_shared<Range>(opt_ents));
765 CHKERR DMSetUp(subset_dm_bdy);
766
767 auto subset_dm_domain = createDM(mField.get_comm(), "DMMOFEM");
768 CHKERR DMMoFEMSetSquareProblem(subset_dm_domain, PETSC_TRUE);
769 CHKERR DMMoFEMCreateSubDM(subset_dm_domain, adjointDM, "SUBSET_DOMAIN");
770 CHKERR DMMoFEMAddElement(subset_dm_domain, "ADJOINT_DOMAIN_FE");
771 CHKERR DMMoFEMAddSubFieldRow(subset_dm_domain, "ADJOINT_FIELD");
772 CHKERR DMMoFEMAddSubFieldCol(subset_dm_domain, "ADJOINT_FIELD");
773 CHKERR DMSetUp(subset_dm_domain);
774
775 // remove dofs on boundary of the domain
776 auto remove_dofs = [&]() {
778
779 std::array<Range, 3> remove_dim_ents;
780 remove_dim_ents[0] =
781 get_range_from_block(mField, "OPT_REMOVE_X", SPACE_DIM - 1);
782 remove_dim_ents[1] =
783 get_range_from_block(mField, "OPT_REMOVE_Y", SPACE_DIM - 1);
784 remove_dim_ents[2] =
785 get_range_from_block(mField, "OPT_REMOVE_Z", SPACE_DIM - 1);
786
787 for (int d = 0; d != 3; ++d) {
788 MOFEM_LOG("WORLD", Sev::inform)
789 << "Removing topology modes on block OPT_REMOVE_" << (char)('X' + d)
790 << " with " << remove_dim_ents[d].size() << " entities";
791 }
792
793 Range body_ents;
794 CHKERR mField.get_moab().get_entities_by_dimension(0, SPACE_DIM, body_ents,
795 true);
796 auto skin = moab::Skinner(&mField.get_moab());
797 Range boundary_ents;
798 CHKERR skin.find_skin(0, body_ents, false, boundary_ents);
799 for (int d = 0; d != 3; ++d) {
800 boundary_ents = subtract(boundary_ents, remove_dim_ents[d]);
801 }
802 ParallelComm *pcomm =
803 ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
804 CHKERR pcomm->filter_pstatus(boundary_ents,
805 PSTATUS_SHARED | PSTATUS_MULTISHARED,
806 PSTATUS_NOT, -1, nullptr);
807 for (auto d = SPACE_DIM - 2; d >= 0; --d) {
808 if (d >= 0) {
809 Range ents;
810 CHKERR mField.get_moab().get_adjacencies(boundary_ents, d, false, ents,
811 moab::Interface::UNION);
812 boundary_ents.merge(ents);
813 } else {
814 Range verts;
815 CHKERR mField.get_moab().get_connectivity(boundary_ents, verts);
816 boundary_ents.merge(verts);
817 }
818 CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
819 boundary_ents);
820 }
821 boundary_ents.merge(opt_ents);
822 CHKERR mField.getInterface<ProblemsManager>()->removeDofsOnEntities(
823 "SUBSET_DOMAIN", "ADJOINT_FIELD", boundary_ents);
824 for (int d = 0; d != 3; ++d) {
825 CHKERR mField.getInterface<ProblemsManager>()->removeDofsOnEntities(
826 "SUBSET_DOMAIN", "ADJOINT_FIELD", remove_dim_ents[d], d, d);
827 }
828
829 // #ifndef NDEBUG
830 if (mField.get_comm_rank() == 0) {
831 CHKERR save_range(mField.get_moab(), "topoMode_boundary_ents.vtk",
832 boundary_ents);
833 }
834 // #endif
835
837 };
838
839 CHKERR remove_dofs();
840
841 auto get_lhs_fe = [&]() {
842 auto fe_lhs = boost::make_shared<BoundaryEle>(mField);
843 fe_lhs->getRuleHook = [](int, int, int p_data) {
844 return 2 * p_data + p_data - 1;
845 };
846 auto &pip = fe_lhs->getOpPtrVector();
848 "GEOMETRY");
851 pip.push_back(new OpMass("ADJOINT_FIELD", "ADJOINT_FIELD",
852 [](double, double, double) { return 1.; }));
853 return fe_lhs;
854 };
855
856 auto get_rhs_fe = [&]() {
857 auto fe_rhs = boost::make_shared<BoundaryEle>(mField);
858 fe_rhs->getRuleHook = [](int, int, int p_data) {
859 return 2 * p_data + p_data - 1;
860 };
861 auto &pip = fe_rhs->getOpPtrVector();
863 "GEOMETRY");
864
865 return fe_rhs;
866 };
867
868 auto block_name = "OPTIMISE";
869 auto mesh_mng = mField.getInterface<MeshsetsManager>();
870 auto bcs = mesh_mng->getCubitMeshsetPtr(
871
872 std::regex((boost::format("%s(.*)") % block_name).str())
873
874 );
875
876 for (auto &v : modeVecs) {
877 v = createDMVector(subset_dm_bdy);
878 }
879
881 struct OpMode : public OP {
882 OpMode(const std::string name,
883 boost::shared_ptr<ObjectiveFunctionData> python_ptr, int id,
884 std::vector<SmartPetscObj<Vec>> mode_vecs,
885 std::vector<std::array<double, 3>> mode_centroids,
886 std::vector<std::array<double, 6>> mode_bboxes, int block_counter,
887 int mode_counter, boost::shared_ptr<Range> range = nullptr)
888 : OP(name, name, OP::OPROW, range), pythonPtr(python_ptr), iD(id),
889 modeVecs(mode_vecs), modeCentroids(mode_centroids),
890 modeBboxes(mode_bboxes), blockCounter(block_counter),
891 modeCounter(mode_counter) {}
892
893 MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
895
896 if (OP::entsPtr) {
897 if (OP::entsPtr->find(this->getFEEntityHandle()) == OP::entsPtr->end())
899 }
900
901 auto nb_rows = data.getIndices().size();
902 if (!nb_rows) {
904 }
905 auto nb_base_functions = data.getN().size2();
906
908 CHKERR pythonPtr->blockModes(iD, OP::getCoordsAtGaussPts(),
909 modeCentroids[blockCounter],
910 modeBboxes[blockCounter], blockModes);
911
912 auto nb_integration_pts = getGaussPts().size2();
913 if (blockModes.size2() != 3 * nb_integration_pts) {
914 MOFEM_LOG("WORLD", Sev::error)
915 << "Number of modes does not match number of integration points: "
916 << blockModes.size2() << "!=" << 3 * nb_integration_pts;
917 CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY, "modes/integration points");
918 }
919
920 VectorDouble nf(nb_rows);
921
922 int nb_modes = blockModes.size1();
923 for (auto mode = 0; mode != nb_modes; ++mode) {
924 nf.clear();
925 // get mode
926 auto t_mode = getFTensor1FromPtr<3>(&blockModes(mode, 0));
927 // get element volume
928 const double vol = OP::getMeasure();
929 // get integration weights
930 auto t_w = OP::getFTensor0IntegrationWeight();
931 // get base function gradient on rows
932 auto t_base = data.getFTensor0N();
933 // loop over integration points
934 for (int gg = 0; gg != nb_integration_pts; gg++) {
935
936 // take into account Jacobian
937 const double alpha = t_w * vol;
938 // loop over rows base functions
939 auto t_nf = getFTensor1FromPtr<SPACE_DIM>(nf.data().data());
940 int rr = 0;
941 for (; rr != nb_rows / SPACE_DIM; ++rr) {
942 t_nf(i) += alpha * t_base * t_mode(i);
943 ++t_base;
944 ++t_nf;
945 }
946 for (; rr < nb_base_functions; ++rr)
947 ++t_base;
948 ++t_w; // move to another integration weight
949 ++t_mode; // move to another mode
950 }
951 Vec vec = modeVecs[modeCounter + mode];
952 auto size = data.getIndices().size();
953 auto *indices = data.getIndices().data().data();
954 auto *nf_data = nf.data().data();
955 CHKERR VecSetValues(vec, size, indices, nf_data, ADD_VALUES);
956 }
957
959 }
960
961 private:
962 boost::shared_ptr<ObjectiveFunctionData> pythonPtr;
963 MatrixDouble blockModes;
964 std::vector<std::array<double, 3>> modeCentroids;
965 std::vector<std::array<double, 6>> modeBboxes;
966 int iD;
967 std::vector<SmartPetscObj<Vec>> modeVecs;
968 int blockCounter;
969 int modeCounter;
970 };
971
972 auto solve_bdy = [&]() {
974
975 auto fe_lhs = get_lhs_fe();
976 auto fe_rhs = get_rhs_fe();
977 int block_counter = 0;
978 int mode_counter = 0;
979 for (auto &bc : bcs) {
980 auto id = bc->getMeshsetId();
981 Range ents;
982 CHKERR mField.get_moab().get_entities_by_handle(bc->getMeshset(), ents,
983 true);
984 auto range = boost::make_shared<Range>(ents);
985 auto &pip_rhs = fe_rhs->getOpPtrVector();
986 pip_rhs.push_back(new OpMode("ADJOINT_FIELD", pythonPtr, id, modeVecs,
987 modeCentroids, modeBBoxes, block_counter,
988 mode_counter, range));
989 CHKERR DMoFEMLoopFiniteElements(subset_dm_bdy, "ADJOINT_BOUNDARY_FE",
990 fe_rhs);
991 pip_rhs.pop_back();
992 int nb_modes;
993 CHKERR pythonPtr->numberOfModes(id, nb_modes);
994 ++block_counter;
995 mode_counter += nb_modes;
996 MOFEM_LOG("WORLD", Sev::inform)
997 << "Setting mode block block: " << bc->getName()
998 << " with ID: " << bc->getMeshsetId()
999 << " total modes: " << mode_counter;
1000 }
1001
1002 for (auto &v : modeVecs) {
1003 CHKERR VecAssemblyBegin(v);
1004 CHKERR VecAssemblyEnd(v);
1005 CHKERR VecGhostUpdateBegin(v, ADD_VALUES, SCATTER_REVERSE);
1006 CHKERR VecGhostUpdateEnd(v, ADD_VALUES, SCATTER_REVERSE);
1007 }
1008
1009 auto M = createDMMatrix(subset_dm_bdy);
1010 fe_lhs->B = M;
1011 CHKERR DMoFEMLoopFiniteElements(subset_dm_bdy, "ADJOINT_BOUNDARY_FE",
1012 fe_lhs);
1013 CHKERR MatAssemblyBegin(M, MAT_FINAL_ASSEMBLY);
1014 CHKERR MatAssemblyEnd(M, MAT_FINAL_ASSEMBLY);
1015
1016 auto solver = createKSP(mField.get_comm());
1017 CHKERR KSPSetOperators(solver, M, M);
1018 CHKERR KSPSetFromOptions(solver);
1019 CHKERR KSPSetUp(solver);
1020 auto v = createDMVector(subset_dm_bdy);
1021 for (auto &f : modeVecs) {
1022 CHKERR KSPSolve(solver, f, v);
1023 CHKERR VecSwap(f, v);
1024 }
1025
1026 for (auto &v : modeVecs) {
1027 CHKERR VecGhostUpdateBegin(v, INSERT_VALUES, SCATTER_FORWARD);
1028 CHKERR VecGhostUpdateEnd(v, INSERT_VALUES, SCATTER_FORWARD);
1029 }
1030
1032 };
1033
1034 CHKERR solve_bdy();
1035
1036 auto get_elastic_fe_lhs = [&]() {
1037 auto fe = boost::make_shared<DomainEle>(mField);
1038 fe->getRuleHook = [](int, int, int p_data) {
1039 return 2 * p_data + p_data - 1;
1040 };
1041 auto &pip = fe->getOpPtrVector();
1043 "GEOMETRY");
1044 CHKERR HookeOps::opFactoryDomainLhs<SPACE_DIM, A, I, DomainEleOp>(
1045 mField, pip, "ADJOINT_FIELD", "MAT_ADJOINT", Sev::noisy);
1046 return fe;
1047 };
1048
1049 auto get_elastic_fe_rhs = [&]() {
1050 auto fe = boost::make_shared<DomainEle>(mField);
1051 fe->getRuleHook = [](int, int, int p_data) {
1052 return 2 * p_data + p_data - 1;
1053 };
1054 auto &pip = fe->getOpPtrVector();
1056 "GEOMETRY");
1057 CHKERR HookeOps::opFactoryDomainRhs<SPACE_DIM, A, I, DomainEleOp>(
1058 mField, pip, "ADJOINT_FIELD", "MAT_ADJOINT", Sev::noisy);
1059 return fe;
1060 };
1061
1062 auto adjoint_gradient_postprocess = [&](auto mode) {
1064 auto post_proc_mesh = boost::make_shared<moab::Core>();
1065 auto post_proc_begin =
1066 boost::make_shared<PostProcBrokenMeshInMoabBaseBegin>(mField,
1067 post_proc_mesh);
1068 auto post_proc_end = boost::make_shared<PostProcBrokenMeshInMoabBaseEnd>(
1069 mField, post_proc_mesh);
1070
1071 auto geom_vec = boost::make_shared<MatrixDouble>();
1072
1073 auto post_proc_fe =
1074 boost::make_shared<PostProcEleDomain>(mField, post_proc_mesh);
1076 post_proc_fe->getOpPtrVector(), {H1}, "GEOMETRY");
1077 post_proc_fe->getOpPtrVector().push_back(
1078 new OpCalculateVectorFieldValues<SPACE_DIM>("ADJOINT_FIELD", geom_vec,
1079 modeVecs[mode]));
1080
1082
1083 post_proc_fe->getOpPtrVector().push_back(
1084
1085 new OpPPMap(
1086
1087 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
1088
1089 {},
1090
1091 {{"MODE", geom_vec}},
1092
1093 {},
1094
1095 {}
1096
1097 )
1098
1099 );
1100
1102 post_proc_begin->getFEMethod());
1103 CHKERR DMoFEMLoopFiniteElements(adjointDM, "ADJOINT_DOMAIN_FE",
1104 post_proc_fe);
1106 post_proc_begin->getFEMethod());
1107
1108 CHKERR post_proc_end->writeFile("mode_" + std::to_string(mode) + ".h5m");
1109
1111 };
1112
1113 auto solve_domain = [&]() {
1115 auto fe_lhs = get_elastic_fe_lhs();
1116 auto fe_rhs = get_elastic_fe_rhs();
1117 auto v = createDMVector(subset_dm_domain);
1118 auto F = vectorDuplicate(v);
1119 fe_rhs->f = F;
1120
1121 auto M = createDMMatrix(subset_dm_domain);
1122 fe_lhs->B = M;
1123 CHKERR DMoFEMLoopFiniteElements(subset_dm_domain, "ADJOINT_DOMAIN_FE",
1124 fe_lhs);
1125 CHKERR MatAssemblyBegin(M, MAT_FINAL_ASSEMBLY);
1126 CHKERR MatAssemblyEnd(M, MAT_FINAL_ASSEMBLY);
1127
1128 auto solver = createKSP(mField.get_comm());
1129 CHKERR KSPSetOperators(solver, M, M);
1130 CHKERR KSPSetFromOptions(solver);
1131 CHKERR KSPSetUp(solver);
1132
1133 int mode_counter = 0;
1134 for (auto &f : modeVecs) {
1135 CHKERR mField.getInterface<FieldBlas>()->setField(0, "ADJOINT_FIELD");
1136 CHKERR DMoFEMMeshToLocalVector(subset_dm_bdy, f, INSERT_VALUES,
1137 SCATTER_REVERSE);
1138 CHKERR VecZeroEntries(F);
1139 CHKERR DMoFEMLoopFiniteElements(subset_dm_domain, "ADJOINT_DOMAIN_FE",
1140 fe_rhs);
1141 CHKERR VecAssemblyBegin(F);
1142 CHKERR VecAssemblyEnd(F);
1143 CHKERR VecGhostUpdateBegin(F, ADD_VALUES, SCATTER_REVERSE);
1144 CHKERR VecGhostUpdateEnd(F, ADD_VALUES, SCATTER_REVERSE);
1145 CHKERR KSPSolve(solver, F, v);
1146 CHKERR VecGhostUpdateBegin(v, INSERT_VALUES, SCATTER_FORWARD);
1147 CHKERR VecGhostUpdateEnd(v, INSERT_VALUES, SCATTER_FORWARD);
1148 CHKERR DMoFEMMeshToLocalVector(subset_dm_domain, v, INSERT_VALUES,
1149 SCATTER_REVERSE);
1150 auto m = createDMVector(adjointDM);
1151 CHKERR DMoFEMMeshToLocalVector(adjointDM, m, INSERT_VALUES,
1152 SCATTER_FORWARD);
1153 f = m;
1154 ++mode_counter;
1155 }
1157 };
1158
1159 CHKERR solve_domain();
1160
1161 for (int i = 0; i < modeVecs.size(); ++i) {
1162 CHKERR adjoint_gradient_postprocess(i);
1163 }
1164
1166}
1167//! [Adjoint modes]
1168
1169//! [Push operators to pipeline]
1172 auto pip = mField.getInterface<PipelineManager>();
1173
1174 //! [Integration rule]
1175 auto integration_rule = [](int, int, int approx_order) {
1176 return 2 * approx_order + 1;
1177 };
1178
1179 auto integration_rule_bc = [](int, int, int approx_order) {
1180 return 2 * approx_order + 1;
1181 };
1182
1184 CHKERR pip->setDomainLhsIntegrationRule(integration_rule);
1185 CHKERR pip->setBoundaryRhsIntegrationRule(integration_rule_bc);
1186 CHKERR pip->setBoundaryLhsIntegrationRule(integration_rule_bc);
1187 //! [Integration rule]
1188
1190 pip->getOpDomainLhsPipeline(), {H1}, "GEOMETRY");
1192 pip->getOpDomainRhsPipeline(), {H1}, "GEOMETRY");
1194 pip->getOpBoundaryRhsPipeline(), {NOSPACE}, "GEOMETRY");
1196 pip->getOpBoundaryLhsPipeline(), {NOSPACE}, "GEOMETRY");
1197
1198 //! [Push domain stiffness matrix to pipeline]
1199 // Add LHS operator for elasticity (stiffness matrix)
1200 CHKERR HookeOps::opFactoryDomainLhs<SPACE_DIM, A, I, DomainEleOp>(
1201 mField, pip->getOpDomainLhsPipeline(), "U", "MAT_ELASTIC", Sev::noisy);
1202 //! [Push domain stiffness matrix to pipeline]
1203
1204 // Add RHS operator for internal forces
1205 CHKERR HookeOps::opFactoryDomainRhs<SPACE_DIM, A, I, DomainEleOp>(
1206 mField, pip->getOpDomainRhsPipeline(), "U", "MAT_ELASTIC", Sev::noisy);
1207
1208 //! [Push Body forces]
1210 pip->getOpDomainRhsPipeline(), mField, "U", Sev::inform);
1211 //! [Push Body forces]
1212
1213 //! [Push natural boundary conditions]
1214 // Add force boundary condition
1216 pip->getOpBoundaryRhsPipeline(), mField, "U", -1, Sev::inform);
1217 // Add case for mix type of BCs
1219 pip->getOpBoundaryLhsPipeline(), mField, "U", Sev::noisy);
1220 //! [Push natural boundary conditions]
1222}
1223//! [Push operators to pipeline]
1224
1225//! [Solve]
1228 auto simple = mField.getInterface<Simple>();
1229 auto dm = simple->getDM();
1230 auto f = createDMVector(dm);
1231 auto d = vectorDuplicate(f);
1232 CHKERR VecZeroEntries(d);
1233 CHKERR DMoFEMMeshToLocalVector(dm, d, INSERT_VALUES, SCATTER_REVERSE);
1234
1235 auto set_essential_bc = [&]() {
1237 // This is low level pushing finite elements (pipelines) to solver
1238
1239 auto ksp_ctx_ptr = getDMKspCtx(dm);
1240 auto pre_proc_rhs = boost::make_shared<FEMethod>();
1241 auto post_proc_rhs = boost::make_shared<FEMethod>();
1242 auto post_proc_lhs = boost::make_shared<FEMethod>();
1243
1244 auto get_pre_proc_hook = [&]() {
1246 {});
1247 };
1248 pre_proc_rhs->preProcessHook = get_pre_proc_hook();
1249
1250 auto get_post_proc_hook_rhs = [this, post_proc_rhs]() {
1252
1254 post_proc_rhs, 1.)();
1256 };
1257
1258 auto get_post_proc_hook_lhs = [this, post_proc_lhs]() {
1260
1262 post_proc_lhs, 1.)();
1264 };
1265
1266 post_proc_rhs->postProcessHook = get_post_proc_hook_rhs;
1267 post_proc_lhs->postProcessHook = get_post_proc_hook_lhs;
1268
1269 ksp_ctx_ptr->getPreProcComputeRhs().push_front(pre_proc_rhs);
1270 ksp_ctx_ptr->getPostProcComputeRhs().push_back(post_proc_rhs);
1271 ksp_ctx_ptr->getPostProcSetOperators().push_back(post_proc_lhs);
1273 };
1274
1275 auto setup_and_solve = [&](auto solver) {
1277 BOOST_LOG_SCOPED_THREAD_ATTR("Timeline", attrs::timer());
1278 MOFEM_LOG("TIMER", Sev::noisy) << "KSPSetUp";
1279 CHKERR KSPSetUp(solver);
1280 MOFEM_LOG("TIMER", Sev::noisy) << "KSPSetUp <= Done";
1281 MOFEM_LOG("TIMER", Sev::noisy) << "KSPSolve";
1282 CHKERR KSPSolve(solver, f, d);
1283 MOFEM_LOG("TIMER", Sev::noisy) << "KSPSolve <= Done";
1285 };
1286
1287 MOFEM_LOG_CHANNEL("TIMER");
1288 MOFEM_LOG_TAG("TIMER", "timer");
1289
1290 CHKERR set_essential_bc();
1291 CHKERR setup_and_solve(kspElastic);
1292
1293 CHKERR VecGhostUpdateBegin(d, INSERT_VALUES, SCATTER_FORWARD);
1294 CHKERR VecGhostUpdateEnd(d, INSERT_VALUES, SCATTER_FORWARD);
1295 CHKERR DMoFEMMeshToLocalVector(dm, d, INSERT_VALUES, SCATTER_REVERSE);
1296
1297 auto evaluate_field_at_the_point = [&]() {
1299
1300 int coords_dim = 3;
1301 std::array<double, 3> field_eval_coords{0.0, 0.0, 0.0};
1302 PetscBool do_eval_field = PETSC_FALSE;
1303 CHKERR PetscOptionsGetRealArray(NULL, NULL, "-field_eval_coords",
1304 field_eval_coords.data(), &coords_dim,
1305 &do_eval_field);
1306
1307 if (do_eval_field) {
1308
1309 vectorFieldPtr = boost::make_shared<MatrixDouble>();
1310 auto field_eval_data =
1311 mField.getInterface<FieldEvaluatorInterface>()->getData<DomainEle>();
1312
1314 ->buildTree<SPACE_DIM>(field_eval_data, simple->getDomainFEName());
1315
1316 field_eval_data->setEvalPoints(field_eval_coords.data(), 1);
1317 auto no_rule = [](int, int, int) { return -1; };
1318 auto field_eval_fe_ptr = field_eval_data->feMethodPtr.lock();
1319 field_eval_fe_ptr->getRuleHook = no_rule;
1320
1321 field_eval_fe_ptr->getOpPtrVector().push_back(
1323
1325 ->evalFEAtThePoint<SPACE_DIM>(
1326 field_eval_coords.data(), 1e-12, simple->getProblemName(),
1327 simple->getDomainFEName(), field_eval_data,
1329 QUIET);
1330
1331 if (vectorFieldPtr->size1()) {
1332 auto t_disp = getFTensor1FromMat<SPACE_DIM>(*vectorFieldPtr);
1333 if constexpr (SPACE_DIM == 2)
1334 MOFEM_LOG("FieldEvaluator", Sev::inform)
1335 << "U_X: " << t_disp(0) << " U_Y: " << t_disp(1);
1336 else
1337 MOFEM_LOG("FieldEvaluator", Sev::inform)
1338 << "U_X: " << t_disp(0) << " U_Y: " << t_disp(1)
1339 << " U_Z: " << t_disp(2);
1340 }
1341
1343 }
1345 };
1346
1347 CHKERR evaluate_field_at_the_point();
1348
1350}
1351//! [Solve]
1352
1353//! [Postprocess results]
1355 SmartPetscObj<Vec> adjoint_vector) {
1357 auto simple = mField.getInterface<Simple>();
1358 auto pip = mField.getInterface<PipelineManager>();
1359 auto det_ptr = boost::make_shared<VectorDouble>();
1360 auto jac_ptr = boost::make_shared<MatrixDouble>();
1361 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
1362 //! [Postprocess clean]
1363 pip->getDomainRhsFE().reset();
1364 pip->getDomainLhsFE().reset();
1365 pip->getBoundaryRhsFE().reset();
1366 pip->getBoundaryLhsFE().reset();
1367 //! [Postprocess clean]
1368
1369 //! [Postprocess initialise]
1370 auto post_proc_mesh = boost::make_shared<moab::Core>();
1371 auto post_proc_begin = boost::make_shared<PostProcBrokenMeshInMoabBaseBegin>(
1372 mField, post_proc_mesh);
1373 auto post_proc_end = boost::make_shared<PostProcBrokenMeshInMoabBaseEnd>(
1374 mField, post_proc_mesh);
1375 //! [Postprocess initialise]
1376
1377 // lambada function to calculate displacement, strain and stress
1378 auto calculate_stress_ops = [&](auto &pip) {
1379 // Add H1 geometry ops
1381
1382 // Use HookeOps commonDataFactory to get matStrainPtr and matCauchyStressPtr
1383 auto common_ptr = HookeOps::commonDataFactory<SPACE_DIM, GAUSS, DomainEle>(
1384 mField, pip, "U", "MAT_ELASTIC", Sev::noisy);
1385
1386 // Store U and GEOMETRY values if needed
1387 auto u_ptr = boost::make_shared<MatrixDouble>();
1388 pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_ptr));
1389 auto x_ptr = boost::make_shared<MatrixDouble>();
1390 pip.push_back(
1391 new OpCalculateVectorFieldValues<SPACE_DIM>("GEOMETRY", x_ptr));
1392
1393 using DataMapMat = std::map<std::string, boost::shared_ptr<MatrixDouble>>;
1394
1395 DataMapMat vec_map = {{"U", u_ptr}, {"GEOMETRY", x_ptr}};
1396 DataMapMat sym_map = {{"STRAIN", common_ptr->getMatStrain()},
1397 {"STRESS", common_ptr->getMatCauchyStress()}};
1398
1399 if (adjoint_vector) {
1400 auto adjoint_ptr = boost::make_shared<MatrixDouble>();
1402 "U", adjoint_ptr, adjoint_vector));
1403 vec_map["ADJOINT"] = adjoint_ptr;
1404 }
1405
1406 // Return what you need: displacements, coordinates, strain, stress
1407 return boost::make_tuple(u_ptr, x_ptr, common_ptr->getMatStrain(),
1408 common_ptr->getMatCauchyStress(), vec_map,
1409 sym_map);
1410 };
1411
1412 auto get_tag_id_on_pmesh = [&](bool post_proc_skin) {
1413 int def_val_int = 0;
1414 Tag tag_mat;
1415 CHKERR mField.get_moab().tag_get_handle(
1416 "MAT_ELASTIC", 1, MB_TYPE_INTEGER, tag_mat,
1417 MB_TAG_CREAT | MB_TAG_SPARSE, &def_val_int);
1418 auto meshset_vec_ptr =
1419 mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(
1420 std::regex((boost::format("%s(.*)") % "MAT_ELASTIC").str()));
1421
1422 Range skin_ents;
1423 std::unique_ptr<Skinner> skin_ptr;
1424 if (post_proc_skin) {
1425 skin_ptr = std::make_unique<Skinner>(&mField.get_moab());
1426 auto boundary_meshset = simple->getBoundaryMeshSet();
1427 CHKERR mField.get_moab().get_entities_by_handle(boundary_meshset,
1428 skin_ents, true);
1429 }
1430
1431 for (auto m : meshset_vec_ptr) {
1432 Range ents_3d;
1433 CHKERR mField.get_moab().get_entities_by_handle(m->getMeshset(), ents_3d,
1434 true);
1435 int const id = m->getMeshsetId();
1436 ents_3d = ents_3d.subset_by_dimension(SPACE_DIM);
1437 if (post_proc_skin) {
1438 Range skin_faces;
1439 CHKERR skin_ptr->find_skin(0, ents_3d, false, skin_faces);
1440 ents_3d = intersect(skin_ents, skin_faces);
1441 }
1442 CHKERR mField.get_moab().tag_clear_data(tag_mat, ents_3d, &id);
1443 }
1444
1445 return tag_mat;
1446 };
1447
1448 auto post_proc_domain = [&](auto post_proc_mesh) {
1449 auto post_proc_fe =
1450 boost::make_shared<PostProcEleDomain>(mField, post_proc_mesh);
1452
1453 auto [u_ptr, x_ptr, mat_strain_ptr, mat_stress_ptr, vec_map, sym_map] =
1454 calculate_stress_ops(post_proc_fe->getOpPtrVector());
1455
1456 post_proc_fe->getOpPtrVector().push_back(
1457
1458 new OpPPMap(
1459
1460 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
1461
1462 {},
1463
1464 vec_map,
1465
1466 {},
1467
1468 sym_map
1469
1470 )
1471
1472 );
1473
1474 post_proc_fe->setTagsToTransfer({get_tag_id_on_pmesh(false)});
1475 return post_proc_fe;
1476 };
1477
1478 auto post_proc_boundary = [&](auto post_proc_mesh) {
1479 auto post_proc_fe =
1480 boost::make_shared<PostProcEleBdy>(mField, post_proc_mesh);
1482 post_proc_fe->getOpPtrVector(), {}, "GEOMETRY");
1483 auto op_loop_side =
1484 new OpLoopSide<SideEle>(mField, simple->getDomainFEName(), SPACE_DIM);
1485 // push ops to side element, through op_loop_side operator
1486 auto [u_ptr, x_ptr, mat_strain_ptr, mat_stress_ptr, vec_map, sym_map] =
1487 calculate_stress_ops(op_loop_side->getOpPtrVector());
1488 post_proc_fe->getOpPtrVector().push_back(op_loop_side);
1489 auto mat_traction_ptr = boost::make_shared<MatrixDouble>();
1490 post_proc_fe->getOpPtrVector().push_back(
1491 new ElasticOps::OpCalculateTraction(mat_stress_ptr, mat_traction_ptr));
1492 vec_map["T"] = mat_traction_ptr;
1493
1495
1496 post_proc_fe->getOpPtrVector().push_back(
1497
1498 new OpPPMap(
1499
1500 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
1501
1502 {},
1503
1504 vec_map,
1505
1506 {},
1507
1508 sym_map
1509
1510 )
1511
1512 );
1513
1514 post_proc_fe->setTagsToTransfer({get_tag_id_on_pmesh(true)});
1515 return post_proc_fe;
1516 };
1517
1518 PetscBool post_proc_skin_only = PETSC_FALSE;
1519 if (SPACE_DIM == 3) {
1520 post_proc_skin_only = PETSC_TRUE;
1521 CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-post_proc_skin_only",
1522 &post_proc_skin_only, PETSC_NULLPTR);
1523 }
1524 if (post_proc_skin_only == PETSC_FALSE) {
1525 pip->getDomainRhsFE() = post_proc_domain(post_proc_mesh);
1526 } else {
1527 pip->getBoundaryRhsFE() = post_proc_boundary(post_proc_mesh);
1528 }
1529
1531 post_proc_begin->getFEMethod());
1532 CHKERR pip->loopFiniteElements();
1534 post_proc_end->getFEMethod());
1535
1536 CHKERR post_proc_end->writeFile("out_elastic_" + std::to_string(iter) +
1537 ".h5m");
1539}
1540//! [Postprocess results]
1541
1542//! [calculateGradient]
1543
1545
1546/**
1547 * @brief Adjoint test operator for elastic problems
1548 *
1549 * This operator implements the adjoint (transpose) of the elastic stiffness operator.
1550 * It's used in the adjoint equation K^T * λ = ∂f/∂u to compute sensitivities.
1551 *
1552 * The operator computes: ∫_Ω ∇v : σ dΩ
1553 * where v are test functions and σ is the Cauchy stress tensor.
1554 *
1555 * This is the transpose of the standard elastic operator ∫_Ω ∇u : C : ∇v dΩ
1556 *
1557 * @tparam SPACE_DIM Spatial dimension (2 or 3)
1558 */
1559template <int SPACE_DIM>
1561 : public DomainBaseOp {
1562
1564
1565 OpAdJointTestOp(const std::string field_name,
1566 boost::shared_ptr<HookeOps::CommonData> comm_ptr)
1567 : OP(field_name, field_name, DomainEleOp::OPROW), commPtr(comm_ptr) {}
1568
1569protected:
1570 boost::shared_ptr<HookeOps::CommonData> commPtr;
1572};
1573
1574/**
1575 * @brief Integration implementation for adjoint test operator
1576 *
1577 * Computes the integral: ∫_Ω ∇v : σ dΩ
1578 * where v are test functions and σ is Cauchy stress tensor
1579 *
1580 * This forms the right-hand side of the adjoint equation K^T * λ = ∂f/∂u
1581 *
1582 * @param row_data Test function data for current element
1583 * @return MoFEMErrorCode Success or error code
1584 */
1585template <int SPACE_DIM>
1588 EntitiesFieldData::EntData &row_data) {
1590 // Define tensor indices for Einstein notation
1597
1598 // Get element geometric information
1599 const double vol = OP::getMeasure();
1600 // Get integration weights for Gauss quadrature
1601 auto t_w = OP::getFTensor0IntegrationWeight();
1602 // Get test function gradients ∇v
1603 auto t_row_grad = row_data.getFTensor1DiffN<SPACE_DIM>();
1604 // Get Cauchy stress tensor σ from forward solution
1605 auto t_cauchy_stress =
1606 getFTensor2SymmetricFromMat<SPACE_DIM>(*(commPtr->getMatCauchyStress()));
1607
1608 // Loop over integration points for numerical integration
1609 for (int gg = 0; gg != OP::nbIntegrationPts; gg++) {
1610 // take into account Jacobian
1611 const double alpha = t_w * vol;
1612
1613 // get rhs vector
1614 auto t_nf = OP::template getNf<SPACE_DIM>();
1615 // loop over rows base functions
1616 int rr = 0;
1617 for (; rr != OP::nbRows / SPACE_DIM; rr++) {
1618 // Variation due to change in geometry (diff base)
1619 t_nf(j) += alpha * t_row_grad(i) * t_cauchy_stress(i, j);
1620
1621 ++t_row_grad;
1622 ++t_nf;
1623 }
1624 for (; rr < OP::nbRowBaseFunctions; ++rr) {
1625 ++t_row_grad;
1626 }
1627
1628 ++t_cauchy_stress;
1629 ++t_w; // move to another integration weight
1630 }
1632}
1633
1634template <int SPACE_DIM>
1636 DomainBaseOp> : public DomainBaseOp {
1637
1639
1641 boost::shared_ptr<HookeOps::CommonData> comm_ptr,
1642 boost::shared_ptr<MatrixDouble> jac,
1643 boost::shared_ptr<MatrixDouble> diff_jac,
1644 boost::shared_ptr<VectorDouble> cof_vals)
1645 : OP(field_name, field_name, DomainEleOp::OPROW), commPtr(comm_ptr),
1646 jac(jac), diffJac(diff_jac), cofVals(cof_vals) {}
1647
1648protected:
1649 boost::shared_ptr<HookeOps::CommonData> commPtr;
1650 boost::shared_ptr<MatrixDouble> jac;
1651 boost::shared_ptr<MatrixDouble> diffJac;
1652 boost::shared_ptr<VectorDouble> cofVals;
1654};
1655
1656template <int DIM> inline auto diff_symmetrize(FTensor::Number<DIM>) {
1657
1658 FTensor::Index<'i', DIM> i;
1659 FTensor::Index<'j', DIM> j;
1660 FTensor::Index<'k', DIM> k;
1661 FTensor::Index<'l', DIM> l;
1662
1664
1665 t_diff(i, j, k, l) = 0;
1666 t_diff(0, 0, 0, 0) = 1;
1667 t_diff(1, 1, 1, 1) = 1;
1668
1669 t_diff(1, 0, 1, 0) = 0.5;
1670 t_diff(1, 0, 0, 1) = 0.5;
1671
1672 t_diff(0, 1, 0, 1) = 0.5;
1673 t_diff(0, 1, 1, 0) = 0.5;
1674
1675 if constexpr (DIM == 3) {
1676 t_diff(2, 2, 2, 2) = 1;
1677
1678 t_diff(2, 0, 2, 0) = 0.5;
1679 t_diff(2, 0, 0, 2) = 0.5;
1680 t_diff(0, 2, 0, 2) = 0.5;
1681 t_diff(0, 2, 2, 0) = 0.5;
1682
1683 t_diff(2, 1, 2, 1) = 0.5;
1684 t_diff(2, 1, 1, 2) = 0.5;
1685 t_diff(1, 2, 1, 2) = 0.5;
1686 t_diff(1, 2, 2, 1) = 0.5;
1687 }
1688
1689 return t_diff;
1690};
1691
1692template <int SPACE_DIM>
1703
1704 auto t_diff_symm = diff_symmetrize(FTensor::Number<SPACE_DIM>());
1705
1706 // get element volume
1707 const double vol = OP::getMeasure();
1708 // get integration weights
1709 auto t_w = OP::getFTensor0IntegrationWeight();
1710 // get Jacobian values
1711 auto t_jac = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*(jac));
1712 // get diff Jacobian values
1713 auto t_diff_jac = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*(diffJac));
1714 // get cofactor values
1715 auto t_cof = getFTensor0FromVec(*(cofVals));
1716 // get base function gradient on rows
1717 auto t_row_grad = row_data.getFTensor1DiffN<SPACE_DIM>();
1718 // get fradient of the field
1719 auto t_grad_u =
1720 getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*(commPtr->matGradPtr));
1721 // get field gradient values
1722 auto t_cauchy_stress =
1723 getFTensor2SymmetricFromMat<SPACE_DIM>(*(commPtr->getMatCauchyStress()));
1724 // material stiffness tensor
1725 auto t_D =
1726 getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(*(commPtr->matDPtr));
1727 // loop over integration points
1728 for (int gg = 0; gg != OP::nbIntegrationPts; gg++) {
1729 // take into account Jacobian
1730 const double alpha = t_w * vol;
1731
1732 auto t_det = determinantTensor(t_jac);
1734 CHKERR invertTensor(t_jac, t_det, t_inv_jac);
1735
1736 // Calculate the variation of the gradient due to geometry change
1738 t_diff_inv_jac(i, j) =
1739 -(t_inv_jac(i, I) * t_diff_jac(I, J)) * t_inv_jac(J, j);
1741 t_diff_grad(i, j) = t_grad_u(i, k) * t_diff_inv_jac(k, j);
1742
1743 // Calculate the variation of the strain tensor
1745 t_diff_strain(i, j) = t_diff_symm(i, j, k, l) * t_diff_grad(k, l);
1746
1747 // Calculate the variation of the stress tensor
1749 t_diff_stress(i, j) = t_D(i, j, k, l) * t_diff_strain(k, l);
1750
1751 // get rhs vector
1752 auto t_nf = OP::template getNf<SPACE_DIM>();
1753 // loop over rows base functions
1754 int rr = 0;
1755 for (; rr != OP::nbRows / SPACE_DIM; rr++) {
1756
1758 t_diff_row_grad(k) = t_row_grad(j) * t_diff_inv_jac(j, k);
1759
1760 // Variation due to change in geometry (diff base)
1761 t_nf(j) += alpha * t_diff_row_grad(i) * t_cauchy_stress(i, j);
1762
1763 // Variation due to change in domain (cofactor)
1764 t_nf(j) += (alpha * t_cof) * t_row_grad(i) * t_cauchy_stress(i, j);
1765
1766 // Variation due to change in stress (diff stress)
1767 t_nf(j) += alpha * t_row_grad(i) * t_diff_stress(i, j);
1768
1769 ++t_row_grad;
1770 ++t_nf;
1771 }
1772 for (; rr < OP::nbRowBaseFunctions; ++rr) {
1773 ++t_row_grad;
1774 }
1775
1776 ++t_grad_u;
1777 ++t_cauchy_stress;
1778 ++t_jac;
1779 ++t_diff_jac;
1780 ++t_cof;
1781 ++t_w; // move to another integration weight
1782 }
1784}
1785
1786/**
1787 * @brief Operator for computing objective function and gradients in topology optimization
1788 *
1789 * This operator interfaces with Python-defined objective functions to:
1790 * 1. Evaluate objective function f(u, x) at current design point
1791 * 2. Compute gradients ∂f/∂u and ∂f/∂x for adjoint sensitivity analysis
1792 *
1793 * The objective function typically depends on:
1794 * - u: displacement field from forward elastic solution
1795 * - x: design variables (topology parameters)
1796 * - σ: stress field computed from displacement
1797 *
1798 * This enables flexible objective function definition (compliance, stress
1799 * constraints, volume constraints, etc.) through Python scripting.
1800 */
1803
1804 OpAdJointObjective(boost::shared_ptr<ObjectiveFunctionData> python_ptr,
1805 boost::shared_ptr<HookeOps::CommonData> comm_ptr,
1806 boost::shared_ptr<MatrixDouble> jac_ptr,
1807 boost::shared_ptr<MatrixDouble> diff_jac,
1808 boost::shared_ptr<VectorDouble> cof_vals,
1809 boost::shared_ptr<MatrixDouble> d_grad_ptr,
1810 boost::shared_ptr<MatrixDouble> u_ptr,
1811 boost::shared_ptr<double> glob_objective_ptr,
1812 boost::shared_ptr<double> glob_objective_grad_ptr)
1813 : OP(NOSPACE, OP::OPSPACE), pythonPtr(python_ptr), commPtr(comm_ptr),
1814 jacPtr(jac_ptr), diffJacPtr(diff_jac), cofVals(cof_vals),
1815 dGradPtr(d_grad_ptr), uPtr(u_ptr), globObjectivePtr(glob_objective_ptr),
1816 globObjectiveGradPtr(glob_objective_grad_ptr) {}
1817
1818 /**
1819 * @brief Compute objective function contributions at element level
1820 *
1821 * Evaluates Python objective function with current displacement and stress
1822 * state, and accumulates global objective value and gradients.
1823 */
1824 MoFEMErrorCode doWork(int side, EntityType type,
1827
1828 // Define tensor indices for calculations
1833
1836
1837 constexpr auto symm_size = (SPACE_DIM * (SPACE_DIM + 1)) / 2;
1838
1839 auto t_diff_symm = diff_symmetrize(FTensor::Number<SPACE_DIM>());
1840
1841 auto nb_gauss_pts = getGaussPts().size2();
1842 auto objective_ptr = boost::make_shared<VectorDouble>(nb_gauss_pts);
1843 auto objective_dstress =
1844 boost::make_shared<MatrixDouble>(symm_size, nb_gauss_pts);
1845 auto objective_dstrain =
1846 boost::make_shared<MatrixDouble>(symm_size, nb_gauss_pts);
1847 auto obj_grad = boost::make_shared<MatrixDouble>(SPACE_DIM, nb_gauss_pts);
1848
1849 auto evaluate_python = [&]() {
1851 auto &coords = OP::getCoordsAtGaussPts();
1852 CHKERR pythonPtr->evalObjectiveFunction(
1853 coords, uPtr, commPtr->getMatCauchyStress(), commPtr->getMatStrain(),
1854 objective_ptr);
1855 CHKERR pythonPtr->evalObjectiveGradientStress(
1856 coords, uPtr, commPtr->getMatCauchyStress(), commPtr->getMatStrain(),
1857 objective_dstress);
1858 CHKERR pythonPtr->evalObjectiveGradientStrain(
1859 coords, uPtr, commPtr->getMatCauchyStress(), commPtr->getMatStrain(),
1860 objective_dstrain);
1861
1862 auto t_grad_u =
1863 getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*(commPtr->matGradPtr));
1864 auto t_D =
1865 getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(*(commPtr->matDPtr));
1866 auto t_jac = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*(jacPtr));
1867 auto t_diff_jac = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*(diffJacPtr));
1868 auto t_cof = getFTensor0FromVec(*(cofVals));
1869 auto t_d_grad = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*(dGradPtr));
1870
1871 auto t_obj = getFTensor0FromVec(*objective_ptr);
1872 auto t_obj_dstress =
1873 getFTensor2SymmetricFromMat<SPACE_DIM>(*objective_dstress);
1874 auto t_obj_dstrain =
1875 getFTensor2SymmetricFromMat<SPACE_DIM>(*objective_dstrain);
1876
1877 auto vol = OP::getMeasure();
1878 auto t_w = getFTensor0IntegrationWeight();
1879 for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
1880
1881 auto t_det = determinantTensor(t_jac);
1883 CHKERR invertTensor(t_jac, t_det, t_inv_jac);
1884
1886 t_diff_inv_jac(i, j) =
1887 -(t_inv_jac(i, I) * t_diff_jac(I, J)) * t_inv_jac(J, j);
1889 t_diff_grad(i, j) = t_grad_u(i, k) * t_diff_inv_jac(k, j);
1890
1892 t_d_strain(i, j) = t_diff_symm(i, j, k, l) * (
1893
1894 t_d_grad(k, l)
1895
1896 +
1897
1898 t_diff_grad(k, l)
1899
1900 );
1901
1902 auto alpha = t_w * vol;
1903
1904 (*globObjectivePtr) += alpha * t_obj;
1905 (*globObjectiveGradPtr) +=
1906 alpha *
1907 (
1908
1909 t_obj_dstress(i, j) * (t_D(i, j, k, l) * t_d_strain(k, l))
1910
1911 +
1912
1913 t_obj_dstrain(i, j) * t_d_strain(i, j)
1914
1915 +
1916
1917 t_obj * t_cof
1918
1919 );
1920
1921 ++t_w;
1922 ++t_jac;
1923 ++t_diff_jac;
1924 ++t_cof;
1925
1926 ++t_obj;
1927 ++t_obj_dstress;
1928 ++t_obj_dstrain;
1929
1930 ++t_grad_u;
1931 ++t_d_grad;
1932 }
1934 };
1935
1936 CHKERR evaluate_python();
1937
1939 }
1940
1941private:
1942 boost::shared_ptr<ObjectiveFunctionData> pythonPtr;
1943 boost::shared_ptr<HookeOps::CommonData> commPtr;
1944 boost::shared_ptr<MatrixDouble> jacPtr;
1945 boost::shared_ptr<MatrixDouble> diffJacPtr;
1946 boost::shared_ptr<VectorDouble> cofVals;
1947 boost::shared_ptr<MatrixDouble> dGradPtr;
1948 boost::shared_ptr<MatrixDouble> uPtr;
1949
1950 boost::shared_ptr<double> globObjectivePtr;
1951 boost::shared_ptr<double> globObjectiveGradPtr;
1952};
1953
1954MoFEMErrorCode Example::calculateGradient(PetscReal *objective_function_value,
1955 Vec objective_function_gradient,
1956 Vec adjoint_vector) {
1957 MOFEM_LOG_CHANNEL("WORLD");
1959 auto simple = mField.getInterface<Simple>();
1960
1961 auto ents = get_range_from_block(mField, "OPTIMISE", SPACE_DIM - 1);
1962
1963 auto get_essential_fe = [this]() {
1964 auto post_proc_rhs = boost::make_shared<FEMethod>();
1965 auto get_post_proc_hook_rhs = [this, post_proc_rhs]() {
1967
1969 post_proc_rhs, 0)();
1971 };
1972 post_proc_rhs->postProcessHook = get_post_proc_hook_rhs;
1973 return post_proc_rhs;
1974 };
1975
1976 auto get_fd_direvative_fe = [&]() {
1977 auto fe = boost::make_shared<DomainEle>(mField);
1978 fe->getRuleHook = [](int, int, int p_data) {
1979 return 2 * p_data + p_data - 1;
1980 };
1981 auto &pip = fe->getOpPtrVector();
1983 // Add RHS operators for internal forces
1984 constexpr bool debug = false;
1985 if constexpr (debug) {
1986 auto common_ptr = HookeOps::commonDataFactory<SPACE_DIM, I, DomainEleOp>(
1987 mField, pip, "U", "MAT_ELASTIC", Sev::noisy);
1988 pip.push_back(
1989 new OpAdJointTestOp<SPACE_DIM, I, DomainBaseOp>("U", common_ptr));
1990 } else {
1991 HookeOps::opFactoryDomainRhs<SPACE_DIM, A, I, DomainEleOp>(
1992 mField, pip, "U", "MAT_ELASTIC", Sev::noisy);
1993 }
1994 return fe;
1995 };
1996
1997 auto calulate_fd_residual = [&](auto eps, auto diff_vec, auto fd_vec) {
1999
2000 constexpr bool debug = false;
2001
2002 auto geom_norm = [](MoFEM::Interface &mField) {
2004 auto field_blas = mField.getInterface<FieldBlas>();
2005 double nrm2 = 0.0;
2006 auto norm2_field = [&](const double val) {
2007 nrm2 += val * val;
2008 return val;
2009 };
2010 CHKERR field_blas->fieldLambdaOnValues(norm2_field, "GEOMETRY");
2011 MPI_Allreduce(MPI_IN_PLACE, &nrm2, 1, MPI_DOUBLE, MPI_SUM,
2012 mField.get_comm());
2013 MOFEM_LOG("WORLD", Sev::inform) << "Geometry norm: " << sqrt(nrm2);
2015 };
2016
2017 if constexpr (debug)
2018 CHKERR geom_norm(mField);
2019
2020 auto initial_current_geometry = createDMVector(adjointDM);
2021 CHKERR mField.getInterface<VecManager>()->setOtherLocalGhostVector(
2022 "ADJOINT", "ADJOINT_FIELD", "GEOMETRY", RowColData::ROW,
2023 initial_current_geometry, INSERT_VALUES, SCATTER_FORWARD);
2024 CHKERR VecAssemblyBegin(initial_current_geometry);
2025 CHKERR VecAssemblyEnd(initial_current_geometry);
2026
2027 if constexpr (debug)
2028 CHKERR geom_norm(mField);
2029
2030 auto perturb_geometry = [&](auto eps, auto diff_vec) {
2032 auto current_geometry = vectorDuplicate(initial_current_geometry);
2033 CHKERR VecCopy(initial_current_geometry, current_geometry);
2034 CHKERR VecAXPY(current_geometry, eps, diff_vec);
2035 CHKERR mField.getInterface<VecManager>()->setOtherLocalGhostVector(
2036 "ADJOINT", "ADJOINT_FIELD", "GEOMETRY", RowColData::ROW,
2037 current_geometry, INSERT_VALUES, SCATTER_REVERSE);
2039 };
2040
2041 auto fe = get_fd_direvative_fe();
2042 auto fp = vectorDuplicate(diff_vec);
2043 auto fm = vectorDuplicate(diff_vec);
2044 auto calc_impl = [&](auto f, auto eps) {
2046 CHKERR VecZeroEntries(f);
2047 fe->f = f;
2048 CHKERR perturb_geometry(eps, diff_vec);
2050 simple->getDomainFEName(), fe);
2051 CHKERR VecAssemblyBegin(f);
2052 CHKERR VecAssemblyEnd(f);
2053 CHKERR VecGhostUpdateBegin(f, ADD_VALUES, SCATTER_REVERSE);
2054 CHKERR VecGhostUpdateEnd(f, ADD_VALUES, SCATTER_REVERSE);
2055 auto post_proc_rhs = get_essential_fe();
2056 post_proc_rhs->f = f;
2058 post_proc_rhs.get());
2060 };
2061 CHKERR calc_impl(fp, eps);
2062 CHKERR calc_impl(fm, -eps);
2063 CHKERR VecWAXPY(fd_vec, -1.0, fm, fp);
2064 CHKERR VecScale(fd_vec, 1.0 / (2.0 * eps));
2065
2066 CHKERR mField.getInterface<VecManager>()->setOtherLocalGhostVector(
2067 "ADJOINT", "ADJOINT_FIELD", "GEOMETRY", RowColData::ROW,
2068 initial_current_geometry, INSERT_VALUES, SCATTER_REVERSE);
2069
2070 if constexpr (debug)
2071 CHKERR geom_norm(mField);
2072
2074 };
2075
2076 auto get_direvative_fe = [&](auto diff_vec) {
2077 auto fe_adjoint = boost::make_shared<DomainEle>(mField);
2078 fe_adjoint->getRuleHook = [](int, int, int p_data) {
2079 return 2 * p_data + p_data - 1;
2080 };
2081 auto &pip = fe_adjoint->getOpPtrVector();
2082
2083 auto jac_ptr = boost::make_shared<MatrixDouble>();
2084 auto det_ptr = boost::make_shared<VectorDouble>();
2085 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
2086 auto diff_jac_ptr = boost::make_shared<MatrixDouble>();
2087 auto cof_ptr = boost::make_shared<VectorDouble>();
2088
2089 using OpCoFactor =
2090 AdJoint<DomainEleOp>::Integration<GAUSS>::OpGetCoFactor<SPACE_DIM>;
2091
2093
2095 "GEOMETRY", jac_ptr));
2097 "U", diff_jac_ptr,
2098 diff_vec)); // Note: that vector is stored on displacemnt vector, that
2099 // why is used here
2100
2101 pip.push_back(new OpCoFactor(jac_ptr, diff_jac_ptr, cof_ptr));
2102
2103 // Add RHS operators for internal forces
2104 auto common_ptr = HookeOps::commonDataFactory<SPACE_DIM, I, DomainEleOp>(
2105 mField, pip, "U", "MAT_ELASTIC", Sev::noisy);
2107 "U", common_ptr, jac_ptr, diff_jac_ptr, cof_ptr));
2108
2109 return fe_adjoint;
2110 };
2111
2112 auto get_objective_fe = [&](auto diff_vec, auto grad_vec,
2113 auto glob_objective_ptr,
2114 auto glob_objective_grad_ptr) {
2115 auto fe_adjoint = boost::make_shared<DomainEle>(mField);
2116 fe_adjoint->getRuleHook = [](int, int, int p_data) {
2117 return 2 * p_data + p_data - 1;
2118 };
2119 auto &pip = fe_adjoint->getOpPtrVector();
2121
2122 auto jac_ptr = boost::make_shared<MatrixDouble>();
2123 auto det_ptr = boost::make_shared<VectorDouble>();
2124 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
2125 auto diff_jac_ptr = boost::make_shared<MatrixDouble>();
2126 auto cof_ptr = boost::make_shared<VectorDouble>();
2127 auto d_grad_ptr = boost::make_shared<MatrixDouble>();
2128 auto u_ptr = boost::make_shared<MatrixDouble>();
2129
2130 using OpCoFactor =
2131 AdJoint<DomainEleOp>::Integration<GAUSS>::OpGetCoFactor<SPACE_DIM>;
2133 "GEOMETRY", jac_ptr));
2135 "U", diff_jac_ptr,
2136 diff_vec)); // Note: that vector is stored on displacemnt vector, that
2137 // why is used here
2138 pip.push_back(new OpCoFactor(jac_ptr, diff_jac_ptr, cof_ptr));
2140 "U", d_grad_ptr, grad_vec));
2141 pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_ptr));
2142
2143 auto common_ptr = HookeOps::commonDataFactory<SPACE_DIM, I, DomainEleOp>(
2144 mField, pip, "U", "MAT_ELASTIC", Sev::noisy);
2145 pip.push_back(new OpAdJointObjective(
2146 pythonPtr, common_ptr, jac_ptr, diff_jac_ptr, cof_ptr, d_grad_ptr,
2147 u_ptr, glob_objective_ptr, glob_objective_grad_ptr));
2148
2149 return fe_adjoint;
2150 };
2151
2152 auto dm = simple->getDM();
2153 auto f = createDMVector(dm);
2154 auto d = vectorDuplicate(f);
2155 auto dm_diff_vec = vectorDuplicate(d);
2156
2157 auto adjoint_fe = get_direvative_fe(dm_diff_vec);
2158 auto glob_objective_ptr = boost::make_shared<double>(0.0);
2159 auto glob_objective_grad_ptr = boost::make_shared<double>(0.0);
2160 auto objective_fe = get_objective_fe(dm_diff_vec, d, glob_objective_ptr,
2161 glob_objective_grad_ptr);
2162
2163 auto set_variance_of_geometry = [&](auto mode, auto mod_vec) {
2165 CHKERR DMoFEMMeshToLocalVector(adjointDM, mod_vec, INSERT_VALUES,
2166 SCATTER_REVERSE);
2167 CHKERR mField.getInterface<VecManager>()->setOtherLocalGhostVector(
2168 simple->getProblemName(), "U", "ADJOINT_FIELD", RowColData::ROW,
2169 dm_diff_vec, INSERT_VALUES, SCATTER_FORWARD);
2170 CHKERR VecGhostUpdateBegin(dm_diff_vec, INSERT_VALUES, SCATTER_FORWARD);
2171 CHKERR VecGhostUpdateEnd(dm_diff_vec, INSERT_VALUES, SCATTER_FORWARD);
2173 };
2174
2175 auto calculate_variance_internal_forces = [&](auto mode, auto mod_vec) {
2177 CHKERR VecZeroEntries(f);
2178 CHKERR VecGhostUpdateBegin(f, INSERT_VALUES, SCATTER_FORWARD);
2179 CHKERR VecGhostUpdateEnd(f, INSERT_VALUES, SCATTER_FORWARD);
2180 adjoint_fe->f = f;
2181 CHKERR DMoFEMLoopFiniteElements(dm, simple->getDomainFEName(), adjoint_fe);
2182 CHKERR VecAssemblyBegin(f);
2183 CHKERR VecAssemblyEnd(f);
2184 CHKERR VecGhostUpdateBegin(f, ADD_VALUES, SCATTER_REVERSE);
2185 CHKERR VecGhostUpdateEnd(f, ADD_VALUES, SCATTER_REVERSE);
2186 auto post_proc_rhs = get_essential_fe();
2187 post_proc_rhs->f = f;
2189 post_proc_rhs.get());
2190 CHKERR VecScale(f, -1.0);
2191
2192#ifndef NDEBUG
2193 constexpr bool debug = false;
2194 if constexpr (debug) {
2195 double norm0;
2196 CHKERR VecNorm(f, NORM_2, &norm0);
2197 auto fd_check = vectorDuplicate(f);
2198 double eps = 1e-5;
2199 CHKERR calulate_fd_residual(eps, dm_diff_vec, fd_check);
2200 double nrm;
2201 CHKERR VecAXPY(fd_check, -1.0, f);
2202 CHKERR VecNorm(fd_check, NORM_2, &nrm);
2203 MOFEM_LOG("WORLD", Sev::inform)
2204 << " FD check for internal forces [ " << mode << " ]: " << nrm
2205 << " / " << norm0 << " ( " << (nrm / norm0) << " )";
2206 }
2207#endif
2208
2210 };
2211
2212 auto calulate_variance_of_displacement = [&](auto mode, auto mod_vec) {
2214 CHKERR KSPSolve(kspElastic, f, d);
2215 CHKERR VecGhostUpdateBegin(d, INSERT_VALUES, SCATTER_FORWARD);
2216 CHKERR VecGhostUpdateEnd(d, INSERT_VALUES, SCATTER_FORWARD);
2218 };
2219
2220 auto calculate_variance_of_objective_function = [&](auto mode, auto mod_vec) {
2222 *glob_objective_ptr = 0.0;
2223 *glob_objective_grad_ptr = 0.0;
2224 CHKERR DMoFEMLoopFiniteElements(dm, simple->getDomainFEName(),
2225 objective_fe);
2226 CHKERR VecSetValue(objective_function_gradient, mode,
2227 *glob_objective_grad_ptr, ADD_VALUES);
2228 std::array<double, 2> array = {*glob_objective_ptr,
2229 *glob_objective_grad_ptr};
2230 MPI_Allreduce(MPI_IN_PLACE, array.data(), 2, MPI_DOUBLE, MPI_SUM,
2231 mField.get_comm());
2232 *glob_objective_ptr = array[0];
2233 *glob_objective_grad_ptr = array[1];
2234 (*objective_function_value) += *glob_objective_ptr;
2235 MOFEM_LOG_CHANNEL("WORLD");
2236 MOFEM_LOG("WORLD", Sev::verbose) << " Objective gradient [ " << mode
2237 << " ]: " << *glob_objective_grad_ptr;
2239 };
2240
2241 CHKERR VecZeroEntries(objective_function_gradient);
2242 CHKERR VecZeroEntries(adjoint_vector);
2243 *objective_function_value = 0.0;
2244 int mode = 0;
2245 for (auto mod_vec : modeVecs) {
2246
2247 CHKERR set_variance_of_geometry(mode, mod_vec);
2248 CHKERR calculate_variance_internal_forces(mode, mod_vec);
2249 CHKERR calulate_variance_of_displacement(mode, mod_vec);
2250 CHKERR calculate_variance_of_objective_function(mode, mod_vec);
2251
2252 CHKERR VecAXPY(adjoint_vector, *glob_objective_grad_ptr, dm_diff_vec);
2253 ++mode;
2254 }
2255
2256 (*objective_function_value) /= modeVecs.size();
2257 MOFEM_LOG("WORLD", Sev::verbose)
2258 << "Objective function: " << *glob_objective_ptr;
2259
2260 CHKERR VecAssemblyBegin(objective_function_gradient);
2261 CHKERR VecAssemblyEnd(objective_function_gradient);
2262
2263 CHKERR VecAssemblyBegin(adjoint_vector);
2264 CHKERR VecAssemblyEnd(adjoint_vector);
2265 CHKERR VecGhostUpdateBegin(adjoint_vector, INSERT_VALUES, SCATTER_FORWARD);
2266 CHKERR VecGhostUpdateEnd(adjoint_vector, INSERT_VALUES, SCATTER_FORWARD);
2267
2269}
2270//! [calculateGradient]
2271
2272static char help[] = "...\n\n";
2273
2274/**
2275 * @brief Main function for topology optimization tutorial using adjoint method
2276 *
2277 * This tutorial demonstrates structural topology optimization using:
2278 * - MoFEM finite element library for elasticity analysis
2279 * - Adjoint method for efficient gradient computation
2280 * - TAO optimization library for design optimization
2281 * - Python interface for flexible objective function definition
2282 *
2283 * Workflow:
2284 * 1. Initialize MoFEM, PETSc and Python environments
2285 * 2. Read mesh and setup finite element problem
2286 * 3. Define objective function via Python interface
2287 * 4. Compute topology optimization modes as design variables
2288 * 5. Run gradient-based optimization using adjoint sensitivities
2289 * 6. Post-process optimized design
2290 *
2291 * The adjoint method enables efficient gradient computation with cost
2292 * independent of the number of design variables, making it suitable
2293 * for large-scale topology optimization problems.
2294 *
2295 * Required input files:
2296 * - Mesh file (.h5m format from CUBIT)
2297 * - Parameter file (param_file.petsc)
2298 * - Objective function (objective_function.py)
2299 *
2300 * @param argc Command line argument count
2301 * @param argv Command line argument values
2302 * @return int Exit code (0 for success)
2303 */
2304int main(int argc, char *argv[]) {
2305
2306 // Initialize Python environment for objective function interface
2307 Py_Initialize();
2308 np::initialize();
2309
2310 // Initialize MoFEM/PETSc and MOAB data structures
2311 const char param_file[] = "param_file.petsc";
2312 MoFEM::Core::Initialize(&argc, &argv, param_file, help);
2313
2314 auto core_log = logging::core::get();
2315 core_log->add_sink(
2317
2318 core_log->add_sink(
2319 LogManager::createSink(LogManager::getStrmSync(), "FieldEvaluator"));
2320 LogManager::setLog("FieldEvaluator");
2321 MOFEM_LOG_TAG("FieldEvaluator", "field_eval");
2322
2323 try {
2324
2325 //! [Register MoFEM discrete manager in PETSc]
2326 DMType dm_name = "DMMOFEM";
2327 CHKERR DMRegister_MoFEM(dm_name);
2328 DMType dm_name_mg = "DMMOFEM_MG";
2330 //! [Register MoFEM discrete manager in PETSc
2331
2332 //! [Create MoAB]
2333 moab::Core mb_instance; ///< mesh database
2334 moab::Interface &moab = mb_instance; ///< mesh database interface
2335 //! [Create MoAB]
2336
2337 //! [Create MoFEM]
2338 MoFEM::Core core(moab); ///< finite element database
2339 MoFEM::Interface &m_field = core; ///< finite element database interface
2340 //! [Create MoFEM]
2341
2342 //! [Example]
2343 Example ex(m_field);
2344 CHKERR ex.runProblem();
2345 //! [Example]
2346 }
2348
2350
2351 if (Py_FinalizeEx() < 0) {
2352 exit(120);
2353 }
2354}
2355
2356/**
2357 * @brief Implementation of ObjectiveFunctionData interface using Python integration
2358 *
2359 * This class provides a concrete implementation of the ObjectiveFunctionData interface
2360 * that bridges MoFEM C++ data structures with Python-defined objective functions.
2361 * It enables flexible definition of optimization objectives through Python scripting
2362 * while maintaining high-performance computation in the C++ finite element framework.
2363 *
2364 * Key features:
2365 * - Python function evaluation with automatic data conversion
2366 * - Support for objective function and gradient computations
2367 * - Topology optimization mode definition through Python
2368 * - Efficient NumPy array interfacing for large data sets
2369 * - Automatic memory management between C++ and Python
2370 *
2371 * The class handles:
2372 * 1. Loading and executing Python objective function scripts
2373 * 2. Converting MoFEM data structures to NumPy arrays
2374 * 3. Calling Python functions for objective evaluation
2375 * 4. Converting Python results back to MoFEM format
2376 * 5. Managing Python interpreter state and namespace
2377 *
2378 * Example Python interface functions that must be defined:
2379 * - objectiveFunction(coords, u, stress, strain) -> objective_value
2380 * - objectiveGradientStress(coords, u, stress, strain) -> gradient_wrt_stress
2381 * - numberOfModes(block_id) -> number_of_topology_modes
2382 * - blockModes(block_id, coords, centroid, bbox) -> mode_vectors
2383 */
2386 virtual ~ObjectiveFunctionDataImpl() = default;
2387
2388 /// Initialize Python interpreter and load objective function script
2389 MoFEMErrorCode initPython(const std::string py_file);
2390
2391 /**
2392 * @brief Evaluate objective function at current state
2393 *
2394 * Calls Python-defined objective function with current displacement,
2395 * stress, and strain fields. Used during optimization to compute
2396 * the objective value that drives the optimization process.
2397 *
2398 * @param coords Gauss point coordinates
2399 * @param u_ptr Displacement field values
2400 * @param stress_ptr Stress tensor values
2401 * @param strain_ptr Strain tensor values
2402 * @param o_ptr Output objective function values
2403 * @return MoFEMErrorCode Success or error code
2404 */
2407 boost::shared_ptr<MatrixDouble> u_ptr,
2408 boost::shared_ptr<MatrixDouble> stress_ptr,
2409 boost::shared_ptr<MatrixDouble> strain_ptr,
2410 boost::shared_ptr<VectorDouble> o_ptr);
2411
2412 /**
2413 * @brief Compute gradient of objective function with respect to stress
2414 *
2415 * Evaluates ∂f/∂σ where f is objective function and σ is stress tensor.
2416 * This gradient is used in the adjoint method to compute sensitivities
2417 * efficiently. Essential for gradient-based topology optimization.
2418 *
2419 * @param coords Gauss point coordinates
2420 * @param u_ptr Displacement field values
2421 * @param stress_ptr Stress tensor values
2422 * @param strain_ptr Strain tensor values
2423 * @param o_ptr Output gradient values
2424 * @return MoFEMErrorCode Success or error code
2425 */
2428 boost::shared_ptr<MatrixDouble> u_ptr,
2429 boost::shared_ptr<MatrixDouble> stress_ptr,
2430 boost::shared_ptr<MatrixDouble> strain_ptr,
2431 boost::shared_ptr<MatrixDouble> o_ptr);
2432
2433 /**
2434 * @brief Compute gradient of objective function with respect to strain
2435 *
2436 * Evaluates ∂f/∂ε where f is objective function and ε is strain tensor.
2437 * Used when objective function depends directly on strain measures,
2438 * complementing stress-based gradients in adjoint sensitivity analysis.
2439 *
2440 * @param coords Gauss point coordinates
2441 * @param u_ptr Displacement field values
2442 * @param stress_ptr Stress tensor values
2443 * @param strain_ptr Strain tensor values
2444 * @param o_ptr Output gradient values
2445 * @return MoFEMErrorCode Success or error code
2446 */
2449 boost::shared_ptr<MatrixDouble> u_ptr,
2450 boost::shared_ptr<MatrixDouble> stress_ptr,
2451 boost::shared_ptr<MatrixDouble> strain_ptr,
2452 boost::shared_ptr<MatrixDouble> o_ptr);
2453
2454 /**
2455 * @brief Compute gradient of objective function with respect to displacement
2456 *
2457 * Evaluates ∂f/∂u where f is objective function and u is displacement field.
2458 * This provides direct sensitivity of objective to displacement changes,
2459 * used as right-hand side in adjoint equation K^T * λ = ∂f/∂u.
2460 *
2461 * @param coords Gauss point coordinates
2462 * @param u_ptr Displacement field values
2463 * @param stress_ptr Stress tensor values
2464 * @param strain_ptr Strain tensor values
2465 * @param o_ptr Output gradient values
2466 * @return MoFEMErrorCode Success or error code
2467 */
2470 boost::shared_ptr<MatrixDouble> u_ptr,
2471 boost::shared_ptr<MatrixDouble> stress_ptr,
2472 boost::shared_ptr<MatrixDouble> strain_ptr,
2473 boost::shared_ptr<MatrixDouble> o_ptr);
2474
2475 /// Return number of topology optimization modes for given material block
2476 MoFEMErrorCode numberOfModes(int block_id, int &modes);
2477
2478 /**
2479 * @brief Define spatial topology modes for design optimization
2480 *
2481 * Generates basis functions that define how the geometry can be modified
2482 * during topology optimization. These modes serve as design variables
2483 * and define the design space for optimization.
2484 *
2485 * @param block_id Material block identifier
2486 * @param coords Element coordinates
2487 * @param centroid Block centroid coordinates
2488 * @param bbodx Bounding box dimensions [xmin,xmax,ymin,ymax,zmin,zmax]
2489 * @param o_ptr Output mode vectors
2490 * @return MoFEMErrorCode Success or error code
2491 */
2492 MoFEMErrorCode blockModes(int block_id, MatrixDouble &coords,
2493 std::array<double, 3> &centroid,
2494 std::array<double, 6> &bbodx, MatrixDouble &o_ptr);
2495
2496private:
2497 // Python interpreter objects for objective function evaluation
2498 bp::object mainNamespace; ///< Main Python namespace for script execution
2499 bp::object objectiveFunction; ///< Python function: f(coords,u,stress,strain) -> objective
2500 bp::object objectiveGradientStress; ///< Python function: ∂f/∂σ(coords,u,stress,strain) -> gradient
2501 bp::object objectiveGradientStrain; ///< Python function: ∂f/∂ε(coords,u,stress,strain) -> gradient
2502 bp::object objectiveGradientU; ///< Python function: ∂f/∂u(coords,u,stress,strain) -> gradient
2503 bp::object objectNumberOfModes; ///< Python function: numberOfModes(block_id) -> int
2504 bp::object objectBlockModes; ///< Python function: blockModes(block_id,coords,centroid,bbox) -> modes
2505
2506 /**
2507 * @brief Internal implementation for objective function evaluation
2508 *
2509 * Handles low-level Python function call with NumPy array conversion.
2510 * Converts MoFEM matrices to NumPy arrays, calls Python function,
2511 * and handles return value conversion.
2512 *
2513 * @param coords NumPy array of coordinates
2514 * @param u NumPy array of displacements
2515 * @param stress NumPy array of stress tensors
2516 * @param strain NumPy array of strain tensors
2517 * @param o Output NumPy array for objective values
2518 * @return MoFEMErrorCode Success or error code
2519 */
2521
2522 np::ndarray coords, np::ndarray u,
2523
2524 np::ndarray stress, np::ndarray strain, np::ndarray &o
2525
2526 );
2527
2528 /**
2529 * @brief Internal implementation for stress gradient computation
2530 *
2531 * Calls Python function to compute ∂f/∂σ with automatic array conversion.
2532 * Essential for adjoint-based sensitivity analysis in topology optimization.
2533 *
2534 * @param coords NumPy array of coordinates
2535 * @param u NumPy array of displacements
2536 * @param stress NumPy array of stress tensors
2537 * @param strain NumPy array of strain tensors
2538 * @param o Output NumPy array for gradient values
2539 * @return MoFEMErrorCode Success or error code
2540 */
2542
2543 np::ndarray coords, np::ndarray u,
2544
2545 np::ndarray stress, np::ndarray strain, np::ndarray &o
2546
2547 );
2548
2549 /**
2550 * @brief Internal implementation for strain gradient computation
2551 *
2552 * Evaluates ∂f/∂ε through Python interface with NumPy array handling.
2553 * Provides strain-based sensitivities for comprehensive gradient computation.
2554 *
2555 * @param coords NumPy array of coordinates
2556 * @param u NumPy array of displacements
2557 * @param stress NumPy array of stress tensors
2558 * @param strain NumPy array of strain tensors
2559 * @param o Output NumPy array for gradient values
2560 * @return MoFEMErrorCode Success or error code
2561 */
2563
2564 np::ndarray coords, np::ndarray u,
2565
2566 np::ndarray stress, np::ndarray strain, np::ndarray &o
2567
2568 );
2569
2570 /**
2571 * @brief Internal implementation for displacement gradient computation
2572 *
2573 * Computes ∂f/∂u through Python interface for adjoint equation right-hand side.
2574 * This gradient drives the adjoint solution that enables efficient sensitivity
2575 * computation independent of design variable count.
2576 *
2577 * @param coords NumPy array of coordinates
2578 * @param u NumPy array of displacements
2579 * @param stress NumPy array of stress tensors
2580 * @param strain NumPy array of strain tensors
2581 * @param o Output NumPy array for gradient values
2582 * @return MoFEMErrorCode Success or error code
2583 */
2585
2586 np::ndarray coords, np::ndarray u,
2587
2588 np::ndarray stress, np::ndarray strain, np::ndarray &o
2589
2590 );
2591
2592 /**
2593 * @brief Internal implementation for topology mode generation
2594 *
2595 * Calls Python function to generate spatial basis functions for topology
2596 * optimization. These modes define the design space and parametrize
2597 * allowable geometry modifications during optimization.
2598 *
2599 * @param block_id Material block identifier
2600 * @param coords NumPy array of element coordinates
2601 * @param centroid NumPy array of block centroid
2602 * @param bbodx NumPy array of bounding box dimensions
2603 * @param o_ptr Output NumPy array for mode vectors
2604 * @return MoFEMErrorCode Success or error code
2605 */
2607
2608 int block_id, np::ndarray coords, np::ndarray centroid, np::ndarray bbodx,
2609 np::ndarray &o_ptr
2610
2611 );
2612
2613 /**
2614 * @brief Convert std::vector to NumPy array for Python interface
2615 *
2616 * Efficient conversion from MoFEM data structures to NumPy arrays
2617 * for seamless Python function calls without data copying.
2618 *
2619 * @param data Source vector data
2620 * @param rows Number of rows in resulting array
2621 * @param nb_gauss_pts Number of Gauss points (affects array structure)
2622 * @return np::ndarray NumPy array for Python use
2623 */
2624 np::ndarray convertToNumPy(std::vector<double> &data, int rows,
2625 int nb_gauss_pts);
2626
2627 /**
2628 * @brief Convert raw pointer to NumPy array for Python interface
2629 *
2630 * Low-level conversion for direct memory access to create NumPy arrays.
2631 * Provides zero-copy conversion when possible for performance.
2632 *
2633 * @param ptr Raw data pointer
2634 * @param s Array size
2635 * @return np::ndarray NumPy array for Python use
2636 */
2637 np::ndarray convertToNumPy(double *ptr, int s);
2638
2639 /// Convert symmetric tensor storage to full matrix format
2641
2642 /// Convert full matrix to symmetric tensor storage format
2643 void copyToSymmetric(double *ptr, MatrixDouble &s);
2644};
2645
2646/**
2647 * @brief Factory function to create Python-integrated objective function interface
2648 *
2649 * Creates and initializes an ObjectiveFunctionDataImpl instance that bridges
2650 * MoFEM finite element computations with Python-defined objective functions.
2651 * This enables flexible objective function definition for topology optimization
2652 * while maintaining computational efficiency.
2653 *
2654 * The Python file must define specific functions with correct signatures:
2655 * - objectiveFunction(coords, u, stress, strain) -> objective_value
2656 * - objectiveGradientStress(coords, u, stress, strain) -> gradient_array
2657 * - objectiveGradientStrain(coords, u, stress, strain) -> gradient_array
2658 * - objectiveGradientU(coords, u, stress, strain) -> gradient_array
2659 * - numberOfModes(block_id) -> integer
2660 * - blockModes(block_id, coords, centroid, bbox) -> mode_array
2661 *
2662 * @param py_file Path to Python file containing objective function definitions
2663 * @return boost::shared_ptr<ObjectiveFunctionData> Configured objective function interface
2664 * @throws MoFEM exception if Python initialization fails
2665 */
2666boost::shared_ptr<ObjectiveFunctionData>
2668 auto ptr = boost::make_shared<ObjectiveFunctionDataImpl>();
2669 CHK_THROW_MESSAGE(ptr->initPython(py_file), "init python");
2670 return ptr;
2671}
2672
2673/**
2674 * @brief Initialize Python interpreter and load objective function script
2675 *
2676 * This method sets up the Python environment for objective function evaluation
2677 * by loading a Python script that defines the required optimization functions.
2678 * It establishes the bridge between MoFEM's C++ finite element computations
2679 * and user-defined Python objective functions.
2680 *
2681 * The Python script must define the following functions:
2682 * - f(coords, u, stress, strain): Main objective function
2683 * - f_stress(coords, u, stress, strain): Gradient w.r.t. stress ∂f/∂σ
2684 * - f_strain(coords, u, stress, strain): Gradient w.r.t. strain ∂f/∂ε
2685 * - f_u(coords, u, stress, strain): Gradient w.r.t. displacement ∂f/∂u
2686 * - number_of_modes(block_id): Return number of topology modes
2687 * - block_modes(block_id, coords, centroid, bbox): Define topology modes
2688 *
2689 * All functions receive NumPy arrays and must return NumPy arrays of
2690 * appropriate dimensions for the finite element computation.
2691 *
2692 * @param py_file Path to Python script containing objective function definitions
2693 * @return MoFEMErrorCode Success or error code
2694 * @throws MOFEM_OPERATION_UNSUCCESSFUL if Python script has errors
2695 */
2697ObjectiveFunctionDataImpl::initPython(const std::string py_file) {
2699 try {
2700
2701 // Create main Python module and namespace for script execution
2702 auto main_module = bp::import("__main__");
2703 mainNamespace = main_module.attr("__dict__");
2704
2705 // Execute the Python script in the main namespace
2706 bp::exec_file(py_file.c_str(), mainNamespace, mainNamespace);
2707
2708 // Create references to required Python functions for later calling
2709 objectiveFunction = mainNamespace["f"]; // Main objective function
2710 objectiveGradientStress = mainNamespace["f_stress"]; // ∂f/∂σ gradient function
2711 objectiveGradientStrain = mainNamespace["f_strain"]; // ∂f/∂ε gradient function
2712 objectiveGradientU = mainNamespace["f_u"]; // ∂f/∂u gradient function
2713 objectNumberOfModes = mainNamespace["number_of_modes"]; // Topology mode count function
2714 objectBlockModes = mainNamespace["block_modes"]; // Topology mode definition function
2715
2716 } catch (bp::error_already_set const &) {
2717 // Handle Python errors by printing to stderr and throwing MoFEM exception
2718 PyErr_Print();
2720 }
2722}
2723
2725 MatrixDouble f(9, s.size2());
2726 f.clear();
2729 auto t_f = getFTensor2FromMat<3, 3>(f);
2730 auto t_s = getFTensor2SymmetricFromMat<SPACE_DIM>(s);
2731 for (int ii = 0; ii != s.size2(); ++ii) {
2732 t_f(i, j) = t_s(i, j);
2733 ++t_f;
2734 ++t_s;
2735 }
2736 return f;
2737};
2738
2743
2744 ptr + 0 * s.size2(), ptr + 1 * s.size2(),
2745 ptr + 2 * s.size2(), ptr + 3 * s.size2(),
2746 ptr + 4 * s.size2(), ptr + 5 * s.size2(),
2747 ptr + 6 * s.size2(), ptr + 7 * s.size2(),
2748 ptr + 8 * s.size2()
2749
2750 };
2751 auto t_s = getFTensor2SymmetricFromMat<SPACE_DIM>(s);
2752 for (int ii = 0; ii != s.size2(); ++ii) {
2753 t_s(i, j) = (t_f(i, j) || t_f(j, i)) / 2.0;
2754 ++t_f;
2755 ++t_s;
2756 }
2757}
2758
2759/**
2760 * @brief Evaluate objective function at current finite element state
2761 *
2762 * This method bridges MoFEM finite element data with Python-defined objective
2763 * functions for topology optimization. It handles the complete data conversion
2764 * workflow from MoFEM matrices to NumPy arrays, calls the Python objective
2765 * function, and converts results back to MoFEM format.
2766 *
2767 * Process:
2768 * 1. Convert coordinate matrix to NumPy format for Python access
2769 * 2. Convert displacement field data to NumPy arrays
2770 * 3. Convert symmetric stress/strain tensors to full 3x3 matrix format
2771 * 4. Call Python objective function: f(coords, u, stress, strain)
2772 * 5. Extract results and copy back to MoFEM vector format
2773 *
2774 * The objective function typically computes scalar quantities like:
2775 * - Compliance: ∫ u^T * f dΩ (minimize structural deformation)
2776 * - Stress constraints: ∫ ||σ - σ_target||² dΩ (control stress distribution)
2777 * - Volume constraints: ∫ ρ dΩ (material usage limitations)
2778 *
2779 * @param coords Gauss point coordinates for current element
2780 * @param u_ptr Displacement field values at Gauss points
2781 * @param stress_ptr Cauchy stress tensor values (symmetric storage)
2782 * @param strain_ptr Strain tensor values (symmetric storage)
2783 * @param o_ptr Output objective function values at each Gauss point
2784 * @return MoFEMErrorCode Success or error code
2785 */
2787 MatrixDouble &coords, boost::shared_ptr<MatrixDouble> u_ptr,
2788 boost::shared_ptr<MatrixDouble> stress_ptr,
2789 boost::shared_ptr<MatrixDouble> strain_ptr,
2790 boost::shared_ptr<VectorDouble> o_ptr) {
2792 try {
2793
2794 // Convert coordinates to NumPy array for Python function
2795 auto np_coords =
2796 convertToNumPy(coords.data(), coords.size1(), coords.size2());
2797 // Convert displacement field to NumPy array
2798 auto np_u = convertToNumPy(u_ptr->data(), u_ptr->size1(), u_ptr->size2());
2799
2800 // Convert symmetric tensor storage to full matrix format for Python
2801 // MoFEM stores symmetric tensors in Voigt notation, Python expects full 3x3 matrices
2802 auto full_stress = copyToFull(*(stress_ptr));
2803 auto full_strain = copyToFull(*(strain_ptr));
2804
2805 // Create NumPy arrays for stress and strain tensors
2806 auto np_stress = convertToNumPy(full_stress.data(), full_stress.size1(),
2807 full_stress.size2());
2808 auto np_strain = convertToNumPy(full_strain.data(), full_strain.size1(),
2809 full_strain.size2());
2810
2811 // Prepare output array for objective function values
2812 np::ndarray np_output = np::empty(bp::make_tuple(strain_ptr->size2()),
2813 np::dtype::get_builtin<double>());
2814
2815 // Call Python objective function implementation
2816 CHKERR objectiveFunctionImpl(np_coords, np_u, np_stress, np_strain,
2817 np_output);
2818
2819 // Copy Python results back to MoFEM vector format
2820 o_ptr->resize(stress_ptr->size2(), false);
2821 double *val_ptr = reinterpret_cast<double *>(np_output.get_data());
2822 std::copy(val_ptr, val_ptr + strain_ptr->size2(), o_ptr->data().begin());
2823
2824 } catch (bp::error_already_set const &) {
2825 // Handle Python errors with detailed error reporting
2826 PyErr_Print();
2828 }
2830}
2831
2832/**
2833 * @brief Compute gradient of objective function with respect to stress tensor
2834 *
2835 * This method evaluates ∂f/∂σ, the partial derivative of the objective function
2836 * with respect to the Cauchy stress tensor. This gradient is fundamental to the
2837 * adjoint method for topology optimization, as it provides the driving force
2838 * for the adjoint equation solution.
2839 *
2840 * Mathematical context:
2841 * The adjoint method requires ∂f/∂σ to compute sensitivities efficiently.
2842 * For stress-based objectives like von Mises stress constraints:
2843 * ∂f/∂σ = ∂/∂σ[∫(σ_vm - σ_target)² dΩ] = 2(σ_vm - σ_target) * ∂σ_vm/∂σ
2844 *
2845 * Process:
2846 * 1. Convert all field data to NumPy format for Python processing
2847 * 2. Call Python function f_stress(coords, u, stress, strain)
2848 * 3. Python returns full 3x3 gradient matrices for each Gauss point
2849 * 4. Convert back to symmetric tensor storage used by MoFEM
2850 * 5. Store results for use in adjoint equation assembly
2851 *
2852 * The resulting gradients drive the adjoint solution that enables efficient
2853 * computation of design sensitivities independent of design variable count.
2854 *
2855 * @param coords Gauss point coordinates
2856 * @param u_ptr Displacement field values
2857 * @param stress_ptr Current stress tensor values (symmetric storage)
2858 * @param strain_ptr Current strain tensor values (symmetric storage)
2859 * @param o_ptr Output stress gradients ∂f/∂σ (symmetric storage)
2860 * @return MoFEMErrorCode Success or error code
2861 */
2863 MatrixDouble &coords, boost::shared_ptr<MatrixDouble> u_ptr,
2864 boost::shared_ptr<MatrixDouble> stress_ptr,
2865 boost::shared_ptr<MatrixDouble> strain_ptr,
2866 boost::shared_ptr<MatrixDouble> o_ptr) {
2868 try {
2869
2870 // Convert coordinates and displacement field to NumPy format
2871 auto np_coords =
2872 convertToNumPy(coords.data(), coords.size1(), coords.size2());
2873 auto np_u = convertToNumPy(u_ptr->data(), u_ptr->size1(), u_ptr->size2());
2874
2875 // Convert symmetric tensors to full 3x3 format for Python processing
2876 auto full_stress = copyToFull(*(stress_ptr));
2877 auto full_strain = copyToFull(*(strain_ptr));
2878
2879 // Create NumPy arrays for stress and strain tensors
2880 auto np_stress = convertToNumPy(full_stress.data(), full_stress.size1(),
2881 full_stress.size2());
2882 auto np_strain = convertToNumPy(full_strain.data(), full_strain.size1(),
2883 full_strain.size2());
2884
2885 // Prepare output array for stress gradients (full matrix format)
2886 np::ndarray np_output =
2887 np::empty(bp::make_tuple(full_strain.size1(), full_strain.size2()),
2888 np::dtype::get_builtin<double>());
2889
2890 // Call Python implementation for stress gradient computation
2891 CHKERR objectiveGradientStressImpl(np_coords, np_u, np_stress, np_strain,
2892 np_output);
2893
2894 // Prepare output matrix in symmetric format
2895 o_ptr->resize(stress_ptr->size1(), stress_ptr->size2(), false);
2896 double *val_ptr = reinterpret_cast<double *>(np_output.get_data());
2897 // Convert full matrix results back to symmetric tensor storage
2898 copyToSymmetric(val_ptr, *(o_ptr));
2899
2900 } catch (bp::error_already_set const &) {
2901 PyErr_Print();
2903 }
2905}
2906
2907/**
2908 * @brief Compute gradient of objective function with respect to strain tensor
2909 *
2910 * This method evaluates ∂f/∂ε, the partial derivative of the objective function
2911 * with respect to the strain tensor. While many structural objectives depend
2912 * primarily on stress, strain-based gradients are important for certain
2913 * optimization formulations and provide additional sensitivity information.
2914 *
2915 * Mathematical context:
2916 * For strain energy-based objectives: f = ½ε:C:ε
2917 * The gradient is: ∂f/∂ε = C:ε = σ (stress tensor)
2918 *
2919 * For strain-based constraints or objectives like strain concentration:
2920 * ∂f/∂ε = ∂/∂ε[∫(ε_vm - ε_target)² dΩ] = 2(ε_vm - ε_target) * ∂ε_vm/∂ε
2921 *
2922 * Process:
2923 * 1. Convert field data to NumPy format for Python compatibility
2924 * 2. Call Python function f_strain(coords, u, stress, strain)
2925 * 3. Python returns gradient matrices ∂f/∂ε for each Gauss point
2926 * 4. Convert from full 3x3 format back to symmetric storage
2927 * 5. Results used in adjoint sensitivity analysis
2928 *
2929 * This gradient complements stress-based gradients in comprehensive
2930 * sensitivity analysis for topology optimization problems.
2931 *
2932 * @param coords Gauss point coordinates
2933 * @param u_ptr Displacement field values
2934 * @param stress_ptr Current stress tensor values (symmetric storage)
2935 * @param strain_ptr Current strain tensor values (symmetric storage)
2936 * @param o_ptr Output strain gradients ∂f/∂ε (symmetric storage)
2937 * @return MoFEMErrorCode Success or error code
2938 */
2940 MatrixDouble &coords, boost::shared_ptr<MatrixDouble> u_ptr,
2941 boost::shared_ptr<MatrixDouble> stress_ptr,
2942 boost::shared_ptr<MatrixDouble> strain_ptr,
2943 boost::shared_ptr<MatrixDouble> o_ptr) {
2945 try {
2946
2947 // Convert coordinates and displacement data to NumPy format
2948 auto np_coords =
2949 convertToNumPy(coords.data(), coords.size1(), coords.size2());
2950 auto np_u = convertToNumPy(u_ptr->data(), u_ptr->size1(), u_ptr->size2());
2951
2952 // Convert symmetric tensor data to full 3x3 matrices for Python
2953 auto full_stress = copyToFull(*(stress_ptr));
2954 auto full_strain = copyToFull(*(strain_ptr));
2955 auto np_stress = convertToNumPy(full_stress.data(), full_stress.size1(),
2956 full_stress.size2());
2957 auto np_strain = convertToNumPy(full_strain.data(), full_strain.size1(),
2958 full_strain.size2());
2959
2960 // Prepare output array for strain gradients
2961 np::ndarray np_output =
2962 np::empty(bp::make_tuple(full_strain.size1(), full_strain.size2()),
2963 np::dtype::get_builtin<double>());
2964
2965 // Call Python implementation for strain gradient computation
2966 CHKERR objectiveGradientStrainImpl(np_coords, np_u, np_stress, np_strain,
2967 np_output);
2968 o_ptr->resize(strain_ptr->size1(), strain_ptr->size2(), false);
2969 double *val_ptr = reinterpret_cast<double *>(np_output.get_data());
2970 copyToSymmetric(val_ptr, *(o_ptr));
2971
2972 } catch (bp::error_already_set const &) {
2973 PyErr_Print();
2975 }
2977}
2978
2979/**
2980 * @brief Compute gradient of objective function with respect to displacement field
2981 *
2982 * This method evaluates ∂f/∂u, the partial derivative of the objective function
2983 * with respect to the displacement field. This gradient is crucial for the adjoint
2984 * method as it forms the right-hand side of the adjoint equation: K^T * λ = ∂f/∂u
2985 *
2986 * Mathematical context:
2987 * The adjoint method solves: K^T * λ = ∂f/∂u
2988 * where λ are the adjoint variables (Lagrange multipliers)
2989 *
2990 * For compliance minimization: f = ½u^T * K * u
2991 * The gradient is: ∂f/∂u = K * u (applied forces)
2992 *
2993 * For displacement-based constraints: f = ||u - u_target||²
2994 * The gradient is: ∂f/∂u = 2(u - u_target)
2995 *
2996 * Process:
2997 * 1. Convert all field data to NumPy format for Python processing
2998 * 2. Call Python function f_u(coords, u, stress, strain)
2999 * 3. Python returns displacement gradients for each component
3000 * 4. Copy results directly (no tensor conversion needed for vectors)
3001 * 5. Results drive adjoint equation solution for sensitivity analysis
3002 *
3003 * This gradient is fundamental to adjoint-based topology optimization,
3004 * enabling efficient sensitivity computation for any number of design variables.
3005 *
3006 * @param coords Gauss point coordinates
3007 * @param u_ptr Displacement field values
3008 * @param stress_ptr Current stress tensor values
3009 * @param strain_ptr Current strain tensor values
3010 * @param o_ptr Output displacement gradients ∂f/∂u
3011 * @return MoFEMErrorCode Success or error code
3012 */
3014 MatrixDouble &coords, boost::shared_ptr<MatrixDouble> u_ptr,
3015 boost::shared_ptr<MatrixDouble> stress_ptr,
3016 boost::shared_ptr<MatrixDouble> strain_ptr,
3017 boost::shared_ptr<MatrixDouble> o_ptr) {
3019 try {
3020
3021 // Convert coordinates and displacement field to NumPy format
3022 auto np_coords =
3023 convertToNumPy(coords.data(), coords.size1(), coords.size2());
3024 auto np_u = convertToNumPy(u_ptr->data(), u_ptr->size1(), u_ptr->size2());
3025
3026 // Convert stress and strain tensors to full matrix format
3027 auto full_stress = copyToFull(*(stress_ptr));
3028 auto full_strain = copyToFull(*(strain_ptr));
3029 auto np_stress = convertToNumPy(full_stress.data(), full_stress.size1(),
3030 full_stress.size2());
3031 auto np_strain = convertToNumPy(full_strain.data(), full_strain.size1(),
3032 full_strain.size2());
3033
3034 // Prepare output array for displacement gradients (same size as displacement field)
3035 np::ndarray np_output =
3036 np::empty(bp::make_tuple(u_ptr->size1(), u_ptr->size2()),
3037 np::dtype::get_builtin<double>());
3038
3039 // Call Python implementation for displacement gradient computation
3040 // Note: This should call objectiveGradientUImpl, not objectiveGradientStrainImpl
3041 CHKERR objectiveGradientUImpl(np_coords, np_u, np_stress, np_strain,
3042 np_output);
3043
3044 // Copy results directly to output matrix (no tensor conversion needed for vectors)
3045 o_ptr->resize(u_ptr->size1(), u_ptr->size2(), false);
3046 double *val_ptr = reinterpret_cast<double *>(np_output.get_data());
3047 std::copy(val_ptr, val_ptr + u_ptr->size1() * u_ptr->size2(),
3048 o_ptr->data().begin());
3049
3050 } catch (bp::error_already_set const &) {
3051 // Handle Python errors with detailed reporting
3052 PyErr_Print();
3054 }
3056}
3057
3058/**
3059 * @brief Generate spatial topology modes for design optimization
3060 *
3061 * This method defines the design parameterization for topology optimization by
3062 * generating spatial basis functions (modes) that describe how the geometry
3063 * can be modified during optimization. These modes serve as design variables
3064 * and define the feasible design space for the optimization problem.
3065 *
3066 * Mathematical context:
3067 * The geometry modification is parameterized as: x_new = x_original + Σ(αᵢ * φᵢ(x))
3068 * where αᵢ are design variables and φᵢ(x) are spatial mode functions
3069 *
3070 * Common mode types:
3071 * - Radial basis functions: φ(x) = exp(-||x-c||²/σ²) for localized changes
3072 * - Polynomial modes: φ(x) = xⁿyᵐzᵖ for global shape changes
3073 * - Sinusoidal modes: φ(x) = sin(kx)cos(ly) for periodic patterns
3074 * - Principal component modes: Derived from geometric sensitivity analysis
3075 *
3076 * Process:
3077 * 1. Query Python function for number of modes for this material block
3078 * 2. Convert coordinate data and geometric information to NumPy format
3079 * 3. Call Python function block_modes(block_id, coords, centroid, bbox)
3080 * 4. Python returns mode vectors for each coordinate at each mode
3081 * 5. Reshape and store modes for use as design variables in optimization
3082 *
3083 * The modes enable efficient design space exploration and gradient-based
3084 * optimization while maintaining geometric feasibility and smoothness.
3085 *
3086 * @param block_id Material block identifier for mode generation
3087 * @param coords Element coordinates where modes are evaluated
3088 * @param centroid Geometric centroid of the material block [x,y,z]
3089 * @param bbodx Bounding box dimensions [xmin,xmax,ymin,ymax,zmin,zmax]
3090 * @param o_ptr Output matrix: modes × (coordinates × spatial_dimension)
3091 * @return MoFEMErrorCode Success or error code
3092 */
3094 int block_id, MatrixDouble &coords, std::array<double, 3> &centroid,
3095 std::array<double, 6> &bbodx, MatrixDouble &o_ptr) {
3097 try {
3098
3099 // Query Python function for number of topology modes for this block
3100 int nb_modes = bp::extract<int>(objectNumberOfModes(block_id));
3101
3102 // Convert coordinate matrix to NumPy format for Python processing
3103 auto np_coords =
3104 convertToNumPy(coords.data(), coords.size1(), coords.size2());
3105
3106 // Convert geometric information to NumPy arrays
3107 auto np_centroid = convertToNumPy(centroid.data(), 3); // Block centroid [x,y,z]
3108 auto np_bbodx = convertToNumPy(bbodx.data(), 6); // Bounding box [xmin,xmax,ymin,ymax,zmin,zmax]
3109
3110 // Prepare output array: [coordinates × modes × spatial_dimensions]
3111 np::ndarray np_output =
3112 np::empty(bp::make_tuple(coords.size1(), nb_modes, 3),
3113 np::dtype::get_builtin<double>());
3114
3115 // Call Python implementation to generate topology modes
3116 CHKERR blockModesImpl(block_id, np_coords, np_centroid, np_bbodx,
3117 np_output);
3118
3119 // Reshape output matrix for MoFEM format: [modes × (coordinates * spatial_dimensions)]
3120 o_ptr.resize(nb_modes, coords.size1() * coords.size2(), false);
3121 double *val_ptr = reinterpret_cast<double *>(np_output.get_data());
3122 // Copy flattened mode data to output matrix
3123 std::copy(val_ptr, val_ptr + coords.size1() * coords.size2() * nb_modes,
3124 o_ptr.data().begin());
3125
3126 } catch (bp::error_already_set const &) {
3127 // Handle Python errors in mode generation
3128 PyErr_Print();
3130 }
3132}
3133
3135
3136 np::ndarray coords, np::ndarray u,
3137
3138 np::ndarray stress, np::ndarray strain, np::ndarray &o
3139
3140) {
3142 try {
3143
3144 // call python function
3145 o = bp::extract<np::ndarray>(objectiveFunction(coords, u, stress, strain));
3146
3147 } catch (bp::error_already_set const &) {
3148 // print all other errors to stderr
3149 PyErr_Print();
3151 }
3153}
3154
3156
3157 np::ndarray coords, np::ndarray u,
3158
3159 np::ndarray stress, np::ndarray strain, np::ndarray &o
3160
3161) {
3163 try {
3164
3165 // call python function
3166 o = bp::extract<np::ndarray>(
3167 objectiveGradientStress(coords, u, stress, strain));
3168
3169 } catch (bp::error_already_set const &) {
3170 // print all other errors to stderr
3171 PyErr_Print();
3173 }
3175}
3176
3178
3179 np::ndarray coords, np::ndarray u,
3180
3181 np::ndarray stress, np::ndarray strain, np::ndarray &o
3182
3183) {
3185 try {
3186
3187 // call python function
3188 o = bp::extract<np::ndarray>(
3189 objectiveGradientStrain(coords, u, stress, strain));
3190
3191 } catch (bp::error_already_set const &) {
3192 // print all other errors to stderr
3193 PyErr_Print();
3195 }
3197}
3198
3200
3201 np::ndarray coords, np::ndarray u,
3202
3203 np::ndarray stress, np::ndarray strain, np::ndarray &o
3204
3205) {
3207 try {
3208
3209 // call python function
3210 o = bp::extract<np::ndarray>(objectiveGradientU(coords, u, stress, strain));
3211
3212 } catch (bp::error_already_set const &) {
3213 // print all other errors to stderr
3214 PyErr_Print();
3216 }
3218}
3219
3221 int &modes) {
3223 try {
3224
3225 modes = bp::extract<int>(objectNumberOfModes(block_id));
3226
3227 } catch (bp::error_already_set const &) {
3228 // print all other errors to stderr
3229 PyErr_Print();
3231 }
3233}
3234
3236 np::ndarray coords,
3237 np::ndarray centroid,
3238 np::ndarray bbodx,
3239 np::ndarray &o) {
3241 try {
3242 // call python function
3243 o = bp::extract<np::ndarray>(
3244 objectBlockModes(block_id, coords, centroid, bbodx));
3245 } catch (bp::error_already_set const &) {
3246 // print all other errors to stderr
3247 PyErr_Print();
3249 }
3251}
3252
3253/**
3254 * @brief Converts a std::vector<double> to a NumPy ndarray.
3255 *
3256 * This function wraps the given vector data into a NumPy array with the
3257 * specified number of rows and Gauss points. The resulting ndarray shares
3258 * memory with the input vector, so changes to one will affect the other.
3259 *
3260 * @param data Reference to the vector containing double values to be converted.
3261 * @param rows Number of rows in the resulting NumPy array.
3262 * @param nb_gauss_pts Number of Gauss points (columns) in the resulting NumPy
3263 * array.
3264 * @return np::ndarray NumPy array view of the input data.
3265 *
3266 * @note
3267 * - `size` specifies the shape of the resulting ndarray as a tuple (rows,
3268 * nb_gauss_pts).
3269 * - `stride` specifies the step size in bytes to move to the next element in
3270 * memory. Here, it is set to sizeof(double), indicating contiguous storage for
3271 * each element.
3272 */
3273inline np::ndarray
3274ObjectiveFunctionDataImpl::convertToNumPy(std::vector<double> &data, int rows,
3275 int nb_gauss_pts) {
3276 auto dtype = np::dtype::get_builtin<double>();
3277 auto size = bp::make_tuple(rows, nb_gauss_pts);
3278 auto stride = bp::make_tuple(nb_gauss_pts * sizeof(double), sizeof(double));
3279 return (np::from_data(data.data(), dtype, size, stride, bp::object()));
3280}
3281
3282inline np::ndarray ObjectiveFunctionDataImpl::convertToNumPy(double *ptr,
3283 int s) {
3284 auto dtype = np::dtype::get_builtin<double>();
3285 auto size = bp::make_tuple(s);
3286 auto stride = bp::make_tuple(sizeof(double));
3287 return (np::from_data(ptr, dtype, size, stride, bp::object()));
3288}
3289
3291 const std::string block_name, int dim) {
3292 Range r;
3293
3294 auto mesh_mng = m_field.getInterface<MeshsetsManager>();
3295 auto bcs = mesh_mng->getCubitMeshsetPtr(
3296
3297 std::regex((boost::format("%s(.*)") % block_name).str())
3298
3299 );
3300
3301 for (auto bc : bcs) {
3302 Range faces;
3303 CHK_MOAB_THROW(bc->getMeshsetIdEntitiesByDimension(m_field.get_moab(), dim,
3304 faces, true),
3305 "get meshset ents");
3306 r.merge(faces);
3307 }
3308
3309 for (auto dd = dim - 1; dd >= 0; --dd) {
3310 if (dd >= 0) {
3311 Range ents;
3312 CHK_MOAB_THROW(m_field.get_moab().get_adjacencies(r, dd, false, ents,
3313 moab::Interface::UNION),
3314 "get adjs");
3315 r.merge(ents);
3316 } else {
3317 Range verts;
3318 CHK_MOAB_THROW(m_field.get_moab().get_connectivity(r, verts),
3319 "get verts");
3320 r.merge(verts);
3321 }
3323 m_field.getInterface<CommInterface>()->synchroniseEntities(r), "comm");
3324 }
3325
3326 return r;
3327};
3328
3329MoFEMErrorCode save_range(moab::Interface &moab, const std::string name,
3330 const Range r) {
3332 auto out_meshset = get_temp_meshset_ptr(moab);
3333 CHKERR moab.add_entities(*out_meshset, r);
3334 CHKERR moab.write_file(name.c_str(), "VTK", "", out_meshset->get_ptr(), 1);
3336};
Calculate traction for linear problem.
Implementation of elastic spring bc.
Natural boundary condition applying pressure from fluid.
Implementation of Hookes operator Hookes for linear elastic problems in MoFEM.
#define MOFEM_LOG_SYNCHRONISE(comm)
Synchronise "SYNC" channel.
Implementation of natural boundary conditions.
Boundary conditions in domain, i.e. body forces.
#define FTENSOR_INDEX(DIM, I)
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition acoustic.cpp:69
auto diff_symmetrize(FTensor::Number< DIM >)
Definition adjoint.cpp:1656
static char help[]
[calculateGradient]
Definition adjoint.cpp:2272
PetscBool is_plane_strain
Definition adjoint.cpp:42
constexpr int SPACE_DIM
[Define dimension]
Definition adjoint.cpp:28
constexpr double poisson_ratio
Poisson's ratio ν
Definition adjoint.cpp:38
constexpr int BASE_DIM
[Constants and material properties]
Definition adjoint.cpp:25
constexpr double shear_modulus_G
Shear modulus G = E/(2(1+ν))
Definition adjoint.cpp:40
constexpr IntegrationType I
Use Gauss quadrature for integration.
Definition adjoint.cpp:33
constexpr double bulk_modulus_K
Bulk modulus K = E/(3(1-2ν))
Definition adjoint.cpp:39
boost::shared_ptr< ObjectiveFunctionData > create_python_objective_function(std::string py_file)
Factory function to create Python-integrated objective function interface.
Definition adjoint.cpp:2667
constexpr AssemblyType A
[Define dimension]
Definition adjoint.cpp:32
Range get_range_from_block(MoFEM::Interface &m_field, const std::string block_name, int dim)
Definition adjoint.cpp:3290
constexpr double young_modulus
[Material properties for linear elasticity]
Definition adjoint.cpp:37
int main()
constexpr double a
static const double eps
constexpr int SPACE_DIM
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
@ QUIET
@ ROW
#define CATCH_ERRORS
Catch errors.
@ MF_EXIST
FieldApproximationBase
approximation base
Definition definitions.h:58
@ LASTBASE
Definition definitions.h:69
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition definitions.h:60
@ DEMKOWICZ_JACOBI_BASE
Definition definitions.h:66
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
@ H1
continuous field
Definition definitions.h:85
@ NOSPACE
Definition definitions.h:83
#define MYPCOMM_INDEX
default communicator number PCOMM
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define CHK_MOAB_THROW(err, msg)
Check error code of MoAB function and throw MoFEM exception.
@ MOFEM_NOT_FOUND
Definition definitions.h:33
@ MOFEM_OPERATION_UNSUCCESSFUL
Definition definitions.h:34
@ MOFEM_DATA_INCONSISTENCY
Definition definitions.h:31
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
PostProcEleByDim< SPACE_DIM >::PostProcEleDomain PostProcEleDomain
PostProcEleByDim< SPACE_DIM >::PostProcEleBdy PostProcEleBdy
@ F
auto integration_rule
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
Definition DMMoFEM.cpp:1113
PetscErrorCode DMMoFEMCreateSubDM(DM subdm, DM dm, const char problem_name[])
Must be called by user to set Sub DM MoFEM data structures.
Definition DMMoFEM.cpp:215
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
Definition DMMoFEM.cpp:488
PetscErrorCode DMMoFEMSetSquareProblem(DM dm, PetscBool square_problem)
set squared problem
Definition DMMoFEM.cpp:450
PetscErrorCode DMMoFEMCreateMoFEM(DM dm, MoFEM::Interface *m_field_ptr, const char problem_name[], const MoFEM::BitRefLevel bit_level, const MoFEM::BitRefLevel bit_mask=MoFEM::BitRefLevel().set())
Must be called by user to set MoFEM data structures.
Definition DMMoFEM.cpp:114
PetscErrorCode DMoFEMPostProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
Definition DMMoFEM.cpp:546
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition DMMoFEM.cpp:514
PetscErrorCode DMMoFEMAddSubFieldRow(DM dm, const char field_name[])
Definition DMMoFEM.cpp:238
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition DMMoFEM.cpp:43
MoFEMErrorCode DMRegister_MGViaApproxOrders(const char sname[])
Register DM for Multi-Grid via approximation orders.
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition DMMoFEM.cpp:576
auto createDMVector(DM dm)
Get smart vector from DM.
Definition DMMoFEM.hpp:1234
PetscErrorCode DMMoFEMAddSubFieldCol(DM dm, const char field_name[])
Definition DMMoFEM.cpp:280
auto createDMMatrix(DM dm)
Get smart matrix from DM.
Definition DMMoFEM.hpp:1191
PetscErrorCode DMoFEMPreProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
Definition DMMoFEM.cpp:536
virtual MoFEMErrorCode add_ents_to_finite_element_by_dim(const EntityHandle entities, const int dim, const std::string name, const bool recursive=true)=0
add entities to finite element
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
virtual MoFEMErrorCode modify_finite_element_add_field_col(const std::string &fe_name, const std::string name_row)=0
set field col which finite element use
virtual MoFEMErrorCode modify_finite_element_add_field_row(const std::string &fe_name, const std::string name_row)=0
set field row which finite element use
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string name_field)=0
set finite element field data
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
virtual MoFEMErrorCode add_ents_to_field_by_dim(const Range &ents, const int dim, const std::string &name, int verb=DEFAULT_VERBOSITY)=0
Add entities to field meshset.
virtual MoFEMErrorCode set_field_order(const EntityHandle meshset, const EntityType type, const std::string &name, const ApproximationOrder order, int verb=DEFAULT_VERBOSITY)=0
Set order approximation of the entities in the field.
IntegrationType
Form integrator integration types.
AssemblyType
[Storage and set boundary conditions]
@ GAUSS
Gaussian quadrature integration.
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
#define MOFEM_LOG(channel, severity)
Log.
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
virtual MoFEMErrorCode loop_dofs(const Problem *problem_ptr, const std::string &field_name, RowColData rc, DofMethod &method, int lower_rank, int upper_rank, int verb=DEFAULT_VERBOSITY)=0
Make a loop over dofs.
MoFEMErrorCode getCubitMeshsetPtr(const int ms_id, const CubitBCType cubit_bc_type, const CubitMeshSets **cubit_meshset_ptr) const
get cubit meshset
FTensor::Index< 'i', SPACE_DIM > i
const double v
phase velocity of light in medium (cm/ns)
FTensor::Index< 'J', DIM1 > J
Definition level_set.cpp:30
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
const FTensor::Tensor2< T, Dim, Dim > Vec
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
implementation of Data Operators for Forces and Sources
Definition Common.hpp:10
auto createKSP(MPI_Comm comm)
PetscErrorCode DMMoFEMSetDestroyProblem(DM dm, PetscBool destroy_problem)
Definition DMMoFEM.cpp:434
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
static const bool debug
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
PetscErrorCode PetscOptionsGetRealArray(PetscOptions *, const char pre[], const char name[], PetscReal dval[], PetscInt *nmax, PetscBool *set)
auto createVectorMPI(MPI_Comm comm, PetscInt n, PetscInt N)
Create MPI Vector.
auto getDMKspCtx(DM dm)
Get KSP context data structure used by DM.
Definition DMMoFEM.hpp:1248
static MoFEMErrorCode invertTensor(FTensor::Tensor2< T1, DIM, DIM > &t, T2 &det, FTensor::Tensor2< T3, DIM, DIM > &inv_t)
static auto determinantTensor(FTensor::Tensor2< T, DIM, DIM > &t)
Calculate the determinant of a tensor of rank DIM.
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
auto get_temp_meshset_ptr(moab::Interface &moab)
Create smart pointer to temporary meshset.
PetscErrorCode TaoSetObjectiveAndGradient(Tao tao, Vec x, PetscReal *f, Vec g, void *ctx)
Sets the objective function value and gradient for a TAO optimization solver.
Definition TaoCtx.cpp:178
auto createDM(MPI_Comm comm, const std::string dm_type_name)
Creates smart DM object.
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
auto createTao(MPI_Comm comm)
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
#define EXECUTABLE_DIMENSION
Definition plastic.cpp:13
PetscBool do_eval_field
Evaluate field.
Definition plastic.cpp:119
ElementsAndOps< SPACE_DIM >::SideEle SideEle
Definition plastic.cpp:61
constexpr auto field_name
static constexpr int approx_order
OpBaseImpl< PETSC, EdgeEleOp > OpBase
Definition radiation.cpp:29
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
[Only used with Hooke equation (linear material model)]
Definition seepage.cpp:56
constexpr double g
FTensor::Index< 'm', 3 > m
Boundary conditions marker.
Definition elastic.cpp:39
[Define entities]
Definition elastic.cpp:38
[Example]
Definition plastic.cpp:216
MoFEMErrorCode boundaryCondition()
Apply essential boundary conditions.
MoFEMErrorCode assembleSystem()
Setup operators in finite element pipeline.
MoFEMErrorCode readMesh()
Read mesh from file and setup meshsets.
SmartPetscObj< KSP > kspElastic
Linear solver for elastic problem.
Definition adjoint.cpp:218
std::vector< SmartPetscObj< Vec > > modeVecs
Topology mode vectors (design variables)
Definition adjoint.cpp:223
FieldApproximationBase base
Choice of finite element basis functions
Definition plot_base.cpp:68
std::vector< std::array< double, 3 > > modeCentroids
Centroids of optimization blocks
Definition adjoint.cpp:224
MoFEMErrorCode topologyModes()
Compute topology optimization modes.
Definition adjoint.cpp:753
SmartPetscObj< Vec > initialGeometry
Initial geometry field.
Definition adjoint.cpp:226
int fieldOrder
Polynomial order for approximation.
Definition adjoint.cpp:215
Example(MoFEM::Interface &m_field)
Definition adjoint.cpp:186
SmartPetscObj< Mat > M
MoFEMErrorCode runProblem()
Main driver function for the optimization process.
MoFEMErrorCode calculateGradient(PetscReal *objective_function_value, Vec objective_function_gradient, Vec adjoint_vector)
Calculate objective function gradient using adjoint method.
Definition adjoint.cpp:1954
MoFEMErrorCode setupAdJoint()
Setup adjoint fields and finite elements
Definition adjoint.cpp:630
MoFEM::Interface & mField
Reference to MoFEM interface.
Definition plastic.cpp:226
std::vector< std::array< double, 6 > > modeBBoxes
Bounding boxes of optimization blocks.
Definition adjoint.cpp:225
boost::shared_ptr< ObjectiveFunctionData > pythonPtr
Interface to Python objective function.
Definition adjoint.cpp:220
SmartPetscObj< DM > adjointDM
Data manager for adjoint problem.
Definition adjoint.cpp:219
MoFEMErrorCode setupProblem()
Setup fields, approximation spaces and DOFs.
MoFEMErrorCode postprocessElastic(int iter, SmartPetscObj< Vec > adjoint_vector=nullptr)
Post-process and output results.
Definition adjoint.cpp:1354
MoFEMErrorCode solveElastic()
Solve forward elastic problem.
Definition adjoint.cpp:1226
boost::shared_ptr< MatrixDouble > vectorFieldPtr
Field values at evaluation points.
Definition adjoint.cpp:194
Add operators pushing bases from local to physical configuration.
Boundary condition manager for finite element problem setup.
Managing BitRefLevels.
MoFEMErrorCode synchroniseEntities(Range &ent, std::map< int, Range > *received_ents, int verb=DEFAULT_VERBOSITY)
synchronize entity range on processors (collective)
virtual moab::Interface & get_moab()=0
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
virtual MoFEMErrorCode add_field(const std::string name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add field.
virtual MPI_Comm & get_comm() const =0
virtual int get_comm_rank() const =0
Core (interface) class.
Definition Core.hpp:82
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition Core.cpp:72
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition Core.cpp:118
Deprecated interface functions.
Definition of the displacement bc data structure.
Definition BCData.hpp:72
Data on single entity (This is passed as argument to DataOperator::doWork)
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
const VectorInt & getIndices() const
Get global indices of degrees of freedom on entity.
Class (Function) to enforce essential constrains on the left hand side diagonal.
Definition Essential.hpp:33
Class (Function) to enforce essential constrains on the right hand side diagonal.
Definition Essential.hpp:41
Class (Function) to enforce essential constrains.
Definition Essential.hpp:25
Basic algebra on fields.
Definition FieldBlas.hpp:21
MoFEMErrorCode fieldLambdaOnValues(OneFieldFunctionOnValues lambda, const std::string field_name, Range *ents_ptr=nullptr)
field lambda
Definition FieldBlas.cpp:21
Field evaluator interface.
static boost::shared_ptr< SinkType > createSink(boost::shared_ptr< std::ostream > stream_ptr, std::string comm_filter)
Create a sink object.
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
static boost::shared_ptr< std::ostream > getStrmSync()
Get the strm sync object.
Interface for managing meshsets containing materials and boundary conditions.
Assembly methods.
Definition Natural.hpp:65
Get field gradients at integration pts for scalar field rank 0, i.e. vector field.
Specialization for double precision vector field values calculation.
Element used to execute operators on side of the element.
Post post-proc data at points from hash maps.
Template struct for dimension-specific finite element types.
PipelineManager interface.
MoFEMErrorCode setDomainRhsIntegrationRule(RuleHookFun rule)
Set integration rule for domain right-hand side finite element.
Problem manager is used to build and partition problems.
Projection of edge entities with one mid-node on hierarchical basis.
Simple interface for fast problem set-up.
Definition Simple.hpp:27
MoFEMErrorCode getOptions()
get options
Definition Simple.cpp:180
MoFEMErrorCode getDM(DM *dm)
Get DM.
Definition Simple.cpp:800
intrusive_ptr for managing petsc objects
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Vector manager is used to create vectors \mofem_vectors.
Implementation of ObjectiveFunctionData interface using Python integration.
Definition adjoint.cpp:2384
MoFEMErrorCode objectiveGradientStrainImpl(np::ndarray coords, np::ndarray u, np::ndarray stress, np::ndarray strain, np::ndarray &o)
Internal implementation for strain gradient computation.
Definition adjoint.cpp:3177
MoFEMErrorCode numberOfModes(int block_id, int &modes)
Return number of topology optimization modes for given material block.
Definition adjoint.cpp:3220
MoFEMErrorCode blockModesImpl(int block_id, np::ndarray coords, np::ndarray centroid, np::ndarray bbodx, np::ndarray &o_ptr)
Internal implementation for topology mode generation.
Definition adjoint.cpp:3235
bp::object objectNumberOfModes
Python function: numberOfModes(block_id) -> int.
Definition adjoint.cpp:2503
MoFEMErrorCode evalObjectiveGradientStrain(MatrixDouble &coords, boost::shared_ptr< MatrixDouble > u_ptr, boost::shared_ptr< MatrixDouble > stress_ptr, boost::shared_ptr< MatrixDouble > strain_ptr, boost::shared_ptr< MatrixDouble > o_ptr)
Compute gradient of objective function with respect to strain.
Definition adjoint.cpp:2939
bp::object objectiveGradientStress
Python function: ∂f/∂σ(coords,u,stress,strain) -> gradient
Definition adjoint.cpp:2500
ObjectiveFunctionDataImpl()=default
MoFEMErrorCode initPython(const std::string py_file)
Initialize Python interpreter and load objective function script.
Definition adjoint.cpp:2697
MoFEMErrorCode blockModes(int block_id, MatrixDouble &coords, std::array< double, 3 > &centroid, std::array< double, 6 > &bbodx, MatrixDouble &o_ptr)
Define spatial topology modes for design optimization.
Definition adjoint.cpp:3093
bp::object objectiveGradientStrain
Python function: ∂f/∂ε(coords,u,stress,strain) -> gradient.
Definition adjoint.cpp:2501
MoFEMErrorCode evalObjectiveGradientU(MatrixDouble &coords, boost::shared_ptr< MatrixDouble > u_ptr, boost::shared_ptr< MatrixDouble > stress_ptr, boost::shared_ptr< MatrixDouble > strain_ptr, boost::shared_ptr< MatrixDouble > o_ptr)
Compute gradient of objective function with respect to displacement.
Definition adjoint.cpp:3013
virtual ~ObjectiveFunctionDataImpl()=default
MoFEMErrorCode evalObjectiveGradientStress(MatrixDouble &coords, boost::shared_ptr< MatrixDouble > u_ptr, boost::shared_ptr< MatrixDouble > stress_ptr, boost::shared_ptr< MatrixDouble > strain_ptr, boost::shared_ptr< MatrixDouble > o_ptr)
Compute gradient of objective function with respect to stress.
Definition adjoint.cpp:2862
bp::object objectiveFunction
Python function: f(coords,u,stress,strain) -> objective.
Definition adjoint.cpp:2499
MoFEMErrorCode evalObjectiveFunction(MatrixDouble &coords, boost::shared_ptr< MatrixDouble > u_ptr, boost::shared_ptr< MatrixDouble > stress_ptr, boost::shared_ptr< MatrixDouble > strain_ptr, boost::shared_ptr< VectorDouble > o_ptr)
Evaluate objective function at current state.
Definition adjoint.cpp:2786
bp::object objectiveGradientU
Python function: ∂f/∂u(coords,u,stress,strain) -> gradient.
Definition adjoint.cpp:2502
MatrixDouble copyToFull(MatrixDouble &s)
Convert symmetric tensor storage to full matrix format.
Definition adjoint.cpp:2724
np::ndarray convertToNumPy(std::vector< double > &data, int rows, int nb_gauss_pts)
Convert std::vector to NumPy array for Python interface.
Definition adjoint.cpp:3274
MoFEMErrorCode objectiveGradientUImpl(np::ndarray coords, np::ndarray u, np::ndarray stress, np::ndarray strain, np::ndarray &o)
Internal implementation for displacement gradient computation.
Definition adjoint.cpp:3199
MoFEMErrorCode objectiveGradientStressImpl(np::ndarray coords, np::ndarray u, np::ndarray stress, np::ndarray strain, np::ndarray &o)
Internal implementation for stress gradient computation.
Definition adjoint.cpp:3155
void copyToSymmetric(double *ptr, MatrixDouble &s)
Convert full matrix to symmetric tensor storage format
Definition adjoint.cpp:2739
bp::object objectBlockModes
Python function: blockModes(block_id,coords,centroid,bbox) -> modes.
Definition adjoint.cpp:2504
bp::object mainNamespace
Main Python namespace for script execution.
Definition adjoint.cpp:2498
MoFEMErrorCode objectiveFunctionImpl(np::ndarray coords, np::ndarray u, np::ndarray stress, np::ndarray strain, np::ndarray &o)
Internal implementation for objective function evaluation.
Definition adjoint.cpp:3134
Abstract interface for Python-defined objective functions.
Definition adjoint.cpp:124
virtual ~ObjectiveFunctionData()=default
virtual MoFEMErrorCode evalObjectiveGradientStress(MatrixDouble &coords, boost::shared_ptr< MatrixDouble > u_ptr, boost::shared_ptr< MatrixDouble > stress_ptr, boost::shared_ptr< MatrixDouble > strain_ptr, boost::shared_ptr< MatrixDouble > o_ptr)=0
Evaluate gradient of objective function w.r.t. stress: ∂f/∂σ
virtual MoFEMErrorCode evalObjectiveGradientStrain(MatrixDouble &coords, boost::shared_ptr< MatrixDouble > u_ptr, boost::shared_ptr< MatrixDouble > stress_ptr, boost::shared_ptr< MatrixDouble > strain_ptr, boost::shared_ptr< MatrixDouble > o_ptr)=0
Evaluate gradient of objective function w.r.t. strain: ∂f/∂ε
virtual MoFEMErrorCode numberOfModes(int block_id, int &modes)=0
Return number of topology optimization modes for given block.
virtual MoFEMErrorCode blockModes(int block_id, MatrixDouble &coords, std::array< double, 3 > &centroid, std::array< double, 6 > &bbox, MatrixDouble &o_ptr)=0
Define spatial modes for topology optimization.
virtual MoFEMErrorCode evalObjectiveFunction(MatrixDouble &coords, boost::shared_ptr< MatrixDouble > u_ptr, boost::shared_ptr< MatrixDouble > stress_ptr, boost::shared_ptr< MatrixDouble > strain_ptr, boost::shared_ptr< VectorDouble > o_ptr)=0
Evaluate objective function f(coords, u, stress, strain)
OpAdJointGradTimesSymTensor(const std::string field_name, boost::shared_ptr< HookeOps::CommonData > comm_ptr, boost::shared_ptr< MatrixDouble > jac, boost::shared_ptr< MatrixDouble > diff_jac, boost::shared_ptr< VectorDouble > cof_vals)
Definition adjoint.cpp:1640
Forward declaration of operator for gradient times symmetric tensor operations.
Definition adjoint.cpp:90
Operator for computing objective function and gradients in topology optimization.
Definition adjoint.cpp:1801
boost::shared_ptr< MatrixDouble > dGradPtr
Definition adjoint.cpp:1947
boost::shared_ptr< double > globObjectiveGradPtr
Definition adjoint.cpp:1951
boost::shared_ptr< ObjectiveFunctionData > pythonPtr
Definition adjoint.cpp:1942
boost::shared_ptr< double > globObjectivePtr
Definition adjoint.cpp:1950
ForcesAndSourcesCore::UserDataOperator OP
Definition adjoint.cpp:1802
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Compute objective function contributions at element level.
Definition adjoint.cpp:1824
boost::shared_ptr< HookeOps::CommonData > commPtr
Definition adjoint.cpp:1943
boost::shared_ptr< VectorDouble > cofVals
Definition adjoint.cpp:1946
boost::shared_ptr< MatrixDouble > uPtr
Definition adjoint.cpp:1948
OpAdJointObjective(boost::shared_ptr< ObjectiveFunctionData > python_ptr, boost::shared_ptr< HookeOps::CommonData > comm_ptr, boost::shared_ptr< MatrixDouble > jac_ptr, boost::shared_ptr< MatrixDouble > diff_jac, boost::shared_ptr< VectorDouble > cof_vals, boost::shared_ptr< MatrixDouble > d_grad_ptr, boost::shared_ptr< MatrixDouble > u_ptr, boost::shared_ptr< double > glob_objective_ptr, boost::shared_ptr< double > glob_objective_grad_ptr)
Definition adjoint.cpp:1804
boost::shared_ptr< MatrixDouble > diffJacPtr
Definition adjoint.cpp:1945
boost::shared_ptr< MatrixDouble > jacPtr
Definition adjoint.cpp:1944
OpAdJointTestOp(const std::string field_name, boost::shared_ptr< HookeOps::CommonData > comm_ptr)
Definition adjoint.cpp:1565
PipelineManager::ElementsAndOpsByDim< 2 >::FaceSideEle SideEle
PipelineManager::ElementsAndOpsByDim< 3 >::FaceSideEle SideEle
auto save_range