v0.16.3
Loading...
Searching...
No Matches
mofem/tutorials/vec-7_shape_optimisation/adjoint.cpp
/**
* @file adjoint.cpp
* @brief Topology optimization using adjoint method with Python objective functions
* @details This tutorial demonstrates:
* - Adjoint-based sensitivity analysis for structural topology optimization
* - Integration with Python for flexible objective function definition
* - Higher-order finite element analysis with geometry fields
* - Automatic differentiation using MoFEM's adjoint capabilities
* @author Anonymous author(s) committing under MIT license
* @example mofem/tutorials/vec-7_shape_optimisation/adjoint.cpp
*
*/
#include <boost/python.hpp>
#include <boost/python/def.hpp>
#include <boost/python/numpy.hpp>
namespace bp = boost::python;
namespace np = boost::python::numpy;
#include <MoFEM.hpp>
using namespace MoFEM;
//! [Constants and material properties]
constexpr int BASE_DIM = 1; ///< Dimension of the base functions
//! [Define dimension]
constexpr int SPACE_DIM =
EXECUTABLE_DIMENSION; ///< Space dimension of problem (2D or 3D), set at compile time
//! [Define dimension]
constexpr AssemblyType A =
AssemblyType::PETSC; ///< Use PETSc for matrix/vector assembly
constexpr IntegrationType I =
IntegrationType::GAUSS; ///< Use Gauss quadrature for integration
//! [Material properties for linear elasticity]
constexpr double young_modulus = 1; ///< Young's modulus E
constexpr double poisson_ratio = 0.3; ///< Poisson's ratio ν
constexpr double bulk_modulus_K =
young_modulus /
(3 * (1 - 2 * poisson_ratio)); ///< Bulk modulus K = E/(3(1-2ν))
constexpr double shear_modulus_G =
young_modulus / (2 * (1 + poisson_ratio)); ///< Shear modulus G = E/(2(1+ν))
PetscBool is_plane_strain =
PETSC_FALSE; ///< Flag for plane strain vs plane stress in 2D
//! [Constants and material properties]
//! [Define finite element types and operators]
using EntData =
EntitiesFieldData::EntData; ///< Entity data for field operations
using DomainEle = PipelineManager::ElementsAndOpsByDim<
SPACE_DIM>::DomainEle; ///< Domain finite elements
using BoundaryEle = PipelineManager::ElementsAndOpsByDim<
SPACE_DIM>::BoundaryEle; ///< Boundary finite elements
using DomainEleOp = DomainEle::UserDataOperator; ///< Domain element operators
using BoundaryEleOp =
BoundaryEle::UserDataOperator; ///< Boundary element operators
//! [Define finite element types and operators]
//! [Boundary condition types]
struct DomainBCs {}; ///< Domain boundary conditions marker
struct BoundaryBCs {}; ///< Boundary conditions marker
//! [Natural boundary condition operators]
using DomainRhsBCs = NaturalBC<DomainEleOp>::Assembly<A>::LinearForm<
I>; ///< Domain RHS natural BCs
using OpDomainRhsBCs =
DomainRhsBCs::OpFlux<DomainBCs, 1, SPACE_DIM>; ///< Domain flux operator
using BoundaryRhsBCs = NaturalBC<BoundaryEleOp>::Assembly<A>::LinearForm<
I>; ///< Boundary RHS natural BCs
using OpBoundaryRhsBCs =
BoundaryRhsBCs::OpFlux<BoundaryBCs, 1,
SPACE_DIM>; ///< Boundary flux operator
using BoundaryLhsBCs = NaturalBC<BoundaryEleOp>::Assembly<A>::BiLinearForm<
I>; ///< Boundary LHS natural BCs
using OpBoundaryLhsBCs =
BoundaryLhsBCs::OpFlux<BoundaryBCs, 1,
SPACE_DIM>; ///< Boundary LHS flux operator
template <int DIM> struct PostProcEleByDim;
template <> struct PostProcEleByDim<2> {
using PostProcEleDomain = PostProcBrokenMeshInMoabBaseCont<DomainEle>;
using PostProcEleBdy = PostProcBrokenMeshInMoabBaseCont<BoundaryEle>;
using SideEle = PipelineManager::ElementsAndOpsByDim<2>::FaceSideEle;
};
template <> struct PostProcEleByDim<3> {
using PostProcEleDomain = PostProcBrokenMeshInMoabBaseCont<DomainEle>;
using PostProcEleBdy = PostProcBrokenMeshInMoabBaseCont<BoundaryEle>;
using SideEle = PipelineManager::ElementsAndOpsByDim<3>::FaceSideEle;
};
// Here you can see how the template is being used
using PostProcEleDomain = PostProcEleByDim<SPACE_DIM>::PostProcEleDomain;
using SideEle = PostProcEleByDim<SPACE_DIM>::SideEle;
using PostProcEleBdy = PostProcEleByDim<SPACE_DIM>::PostProcEleBdy;
// Forward declaration
template <int SPACE_DIM, IntegrationType I, typename OpBase>
struct OpAdJointTestOp;
/// Forward declaration of operator for gradient times symmetric tensor operations
template <int SPACE_DIM, IntegrationType I, typename OpBase>
struct OpAdJointGradTimesSymTensor;
enum SensitivityMethod { DIRECT, ADJOINT };
SensitivityMethod derivative_type = ADJOINT;
#include <ElasticSpring.hpp>
#include <FluidLevel.hpp>
#include <CalculateTraction.hpp>
#include <NaturalDomainBC.hpp>
#include <NaturalBoundaryBC.hpp>
#include <HookeOps.hpp>
#include <ElasticPostProc.hpp>
#include <ObjectiveFunctionData.hpp>
using namespace ShapeOptimization;
Range get_range_from_block(MoFEM::Interface &m_field,
const std::string block_name, int dim);
MoFEMErrorCode save_range(moab::Interface &moab, const std::string name,
const Range r);
struct Example {
Example(MoFEM::Interface &m_field) : mField(m_field) {}
/// Main driver function for the optimization process
MoFEMErrorCode runProblem();
private:
MoFEM::Interface &mField; ///< Reference to MoFEM interface
// Differentiate the same quadrature-discretised problem that is solved.
static int integrationRule(int, int, int order) { return 2 * order + 1; }
// Problem setup methods
MoFEMErrorCode readMesh(); ///< Read mesh from file and setup meshsets
MoFEMErrorCode
setupProblem(); ///< Setup fields, approximation spaces and DOFs
MoFEMErrorCode setupAdJoint(); ///< Setup adjoint fields and finite elements
MoFEMErrorCode boundaryCondition(); ///< Apply essential boundary conditions
MoFEMErrorCode topologyModes(); ///< Compute topology optimization modes
MoFEMErrorCode
assembleSystem(); ///< Setup operators in finite element pipeline
// Analysis methods
MoFEMErrorCode solveElastic(); ///< Solve forward elastic problem
MoFEMErrorCode
postprocessElastic(int iter,
SmartPetscObj<Vec> adjoint_vector =
nullptr); ///< Post-process and output results
/// Calculate objective function gradient using adjoint method
MoFEMErrorCode calculateGradient(PetscReal *objective_function_value,
Vec objective_function_gradient,
Vec adjoint_vector);
// Problem configuration
FieldApproximationBase base; ///< Choice of finite element basis functions
int fieldOrder = 2; ///< Polynomial order for approximation
// Solver and data management
SmartPetscObj<KSP> kspElastic; ///< Linear solver for elastic problem
SmartPetscObj<DM> adjointDM; ///< Data manager for adjoint problem
boost::shared_ptr<ObjectiveFunctionData>
pythonPtr; ///< Interface to Python objective function
// Topology optimization data
std::vector<SmartPetscObj<Vec>>
modeVecs; ///< Topology mode vectors (design variables)
std::vector<std::array<double, 3>>
modeCentroids; ///< Centroids of optimization blocks
std::vector<std::array<double, 6>>
modeBBoxes; ///< Bounding boxes of optimization blocks
SmartPetscObj<Vec> initialGeometry; ///< Initial geometry field
};
//! [Run topology optimization problem]
/**
* @brief Main driver for topology optimization using adjoint sensitivity analysis
*
* This function orchestrates the complete topology optimization workflow:
* 1. Initialize Python objective function interface
* 2. Setup finite element problems (forward and adjoint)
* 3. Compute initial elastic solution
* 4. Generate topology optimization modes
* 5. Run TAO optimization loop with adjoint-based gradients
*
* The optimization uses TAO (Toolkit for Advanced Optimization) with L-BFGS
* algorithm. At each iteration:
* - Update geometry based on current design variables
* - Solve forward elastic problem
* - Compute objective function and gradients using adjoint method
* - Post-process results
*
* @return MoFEMErrorCode Success or error code
*/
MoFEMErrorCode Example::runProblem() {
MoFEMFunctionBegin;
// Read objective function from Python file
char objective_function_file_name[255] = "objective_function.py";
CHKERR PetscOptionsGetString(
PETSC_NULLPTR, PETSC_NULLPTR, "-objective_function",
objective_function_file_name, 255, PETSC_NULLPTR);
// Verify that the Python objective function file exists
auto file_exists = [](std::string myfile) {
std::ifstream file(myfile.c_str());
if (file) {
return true;
}
return false;
};
if (!file_exists(objective_function_file_name)) {
MOFEM_LOG("WORLD", Sev::error) << "Objective function file NOT found: "
<< objective_function_file_name;
CHK_THROW_MESSAGE(MOFEM_NOT_FOUND, "file NOT found");
}
// Create Python interface for objective function
pythonPtr = create_python_objective_function(objective_function_file_name);
char sensitivity_method_name[32] = "adjoint";
CHKERR PetscOptionsGetString(PETSC_NULLPTR, PETSC_NULLPTR,
"-sensitivity_method", sensitivity_method_name,
sizeof(sensitivity_method_name), PETSC_NULLPTR);
std::string sensitivity_method = sensitivity_method_name;
std::transform(sensitivity_method.begin(), sensitivity_method.end(),
sensitivity_method.begin(),
[](unsigned char c) { return std::tolower(c); });
if (sensitivity_method == "direct") {
derivative_type = DIRECT;
} else if (sensitivity_method == "adjoint") {
derivative_type = ADJOINT;
} else {
SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
"Unknown -sensitivity_method. Use 'direct' or 'adjoint'.");
}
MOFEM_LOG("WORLD", Sev::inform)
<< "Sensitivity method: " << sensitivity_method;
// Setup finite element problems
CHKERR readMesh(); // Read mesh and meshsets
CHKERR setupProblem(); // Setup displacement field and geometry field
CHKERR setupAdJoint(); // Setup adjoint field for sensitivity analysis
CHKERR boundaryCondition(); // Apply essential boundary conditions
CHKERR assembleSystem(); // Setup finite element operators
// Create linear solver for elastic problem
auto pip = mField.getInterface<PipelineManager>();
kspElastic = pip->createKSP();
CHKERR KSPSetFromOptions(kspElastic);
// Solve initial elastic problem
CHKERR solveElastic();
CHKERR postprocessElastic(-1); // Post-process initial solution
auto create_vec_modes = [&](auto block_name) {
MoFEMFunctionBegin;
auto mesh_mng = mField.getInterface<MeshsetsManager>();
auto bcs = mesh_mng->getCubitMeshsetPtr(
std::regex((boost::format("%s(.*)") % block_name).str())
);
int nb_total_modes = 0;
for (auto &bc : bcs) {
auto id = bc->getMeshsetId();
int nb_modes;
CHKERR pythonPtr->numberOfModes(id, nb_modes);
nb_total_modes += nb_modes;
}
MOFEM_LOG("WORLD", Sev::inform)
<< "Total number of modes to apply: " << nb_total_modes;
modeVecs.resize(nb_total_modes);
MoFEMFunctionReturn(0);
};
auto get_modes_bounding_boxes = [&](auto block_name) {
MoFEMFunctionBegin;
auto mesh_mng = mField.getInterface<MeshsetsManager>();
auto bcs = mesh_mng->getCubitMeshsetPtr(
std::regex((boost::format("%s(.*)") % block_name).str())
);
for (auto &bc : bcs) {
auto meshset = bc->getMeshset();
Range ents;
CHKERR mField.get_moab().get_entities_by_handle(meshset, ents, true);
Range verts;
CHKERR mField.get_moab().get_connectivity(ents, verts, false);
std::vector<double> x(verts.size());
std::vector<double> y(verts.size());
std::vector<double> z(verts.size());
CHKERR mField.get_moab().get_coords(verts, x.data(), y.data(), z.data());
std::array<double, 3> centroid = {0, 0, 0};
for (int i = 0; i != verts.size(); ++i) {
centroid[0] += x[i];
centroid[1] += y[i];
centroid[2] += z[i];
}
MPI_Allreduce(MPI_IN_PLACE, centroid.data(), 3, MPI_DOUBLE, MPI_SUM,
mField.get_comm());
int nb_vertex = verts.size();
MPI_Allreduce(MPI_IN_PLACE, &nb_vertex, 1, MPI_INT, MPI_SUM,
mField.get_comm());
if (nb_vertex) {
centroid[0] /= nb_vertex;
centroid[1] /= nb_vertex;
centroid[2] /= nb_vertex;
}
std::array<double, 6> bbox = {centroid[0], centroid[1], centroid[2],
centroid[0], centroid[1], centroid[2]};
for (int i = 0; i != verts.size(); ++i) {
bbox[0] = std::min(bbox[0], x[i]);
bbox[1] = std::min(bbox[1], y[i]);
bbox[2] = std::min(bbox[2], z[i]);
bbox[3] = std::max(bbox[3], x[i]);
bbox[4] = std::max(bbox[4], y[i]);
bbox[5] = std::max(bbox[5], z[i]);
}
MPI_Allreduce(MPI_IN_PLACE, &bbox[0], 3, MPI_DOUBLE, MPI_MIN,
mField.get_comm());
MPI_Allreduce(MPI_IN_PLACE, &bbox[3], 3, MPI_DOUBLE, MPI_MAX,
mField.get_comm());
MOFEM_LOG("WORLD", Sev::inform)
<< "Block: " << bc->getName() << " centroid: " << centroid[0] << " "
<< centroid[1] << " " << centroid[2];
MOFEM_LOG("WORLD", Sev::inform)
<< "Block: " << bc->getName() << " bbox: " << bbox[0] << " "
<< bbox[1] << " " << bbox[2] << " " << bbox[3] << " " << bbox[4]
<< " " << bbox[5];
modeCentroids.push_back(centroid);
modeBBoxes.push_back(bbox);
}
MoFEMFunctionReturn(0);
};
auto eval_objective_and_gradient = [](Tao tao, Vec x, PetscReal *f, Vec g,
void *ctx) -> PetscErrorCode {
MoFEMFunctionBegin;
auto *ex_ptr = (Example *)ctx;
int iter;
CHKERR TaoGetIterationNumber(tao, &iter);
auto set_geometry = [&](Vec x) {
MoFEMFunctionBegin;
VecScatter ctx;
Vec x_local;
CHKERR VecScatterCreateToAll(x, &ctx, &x_local);
// scatter as many times as you need
CHKERR VecScatterBegin(ctx, x, x_local, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecScatterEnd(ctx, x, x_local, INSERT_VALUES, SCATTER_FORWARD);
// destroy scatter context and local vector when no longer needed
CHKERR VecScatterDestroy(&ctx);
auto current_geometry = vectorDuplicate(ex_ptr->initialGeometry);
CHKERR VecCopy(ex_ptr->initialGeometry, current_geometry);
const double *a;
CHKERR VecGetArrayRead(x_local, &a);
const double *coeff = a;
for (auto &mode_vec : ex_ptr->modeVecs) {
MOFEM_LOG("WORLD", Sev::verbose)
<< "Adding mode with coeff: " << *coeff;
CHKERR VecAXPY(current_geometry, (*coeff), mode_vec);
++coeff;
}
CHKERR VecRestoreArrayRead(x_local, &a);
CHKERR VecGhostUpdateBegin(current_geometry, INSERT_VALUES,
SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(current_geometry, INSERT_VALUES,
SCATTER_FORWARD);
CHKERR ex_ptr->mField.getInterface<VecManager>()
->setOtherLocalGhostVector("ADJOINT", "ADJOINT_FIELD", "GEOMETRY",
RowColData::ROW, current_geometry,
INSERT_VALUES, SCATTER_REVERSE);
CHKERR VecDestroy(&x_local);
MoFEMFunctionReturn(0);
};
CHKERR set_geometry(x);
CHKERR KSPReset(ex_ptr->kspElastic);
CHKERR ex_ptr->solveElastic();
auto simple = ex_ptr->mField.getInterface<Simple>();
auto adjoint_vector = createDMVector(simple->getDM());
CHKERR ex_ptr->calculateGradient(f, g, adjoint_vector);
CHKERR ex_ptr->postprocessElastic(iter, adjoint_vector);
MoFEMFunctionReturn(0);
};
// Create topology optimization modes and setup TAO solver
CHKERR create_vec_modes("OPTIMISE");
CHKERR get_modes_bounding_boxes("OPTIMISE");
/**
* Setup TAO (Toolkit for Advanced Optimization) solver for topology optimization
*
* TAO provides various optimization algorithms. Here we use TAOLMVM (Limited Memory
* Variable Metric) which is a quasi-Newton method suitable for unconstrained
* optimization with gradient information.
*
* The optimization variables are coefficients for the topology modes, and
* gradients are computed using the adjoint method for efficiency.
*/
auto tao = createTao(mField.get_comm());
CHKERR TaoSetType(tao, TAOLMVM);
auto rank = mField.get_comm_rank();
auto g = createVectorMPI(mField.get_comm(), (!rank) ? modeVecs.size() : 0,
PETSC_DECIDE);
CHKERR TaoSetObjectiveAndGradient(tao, g, eval_objective_and_gradient,
(void *)this);
// Store initial geometry for reference during optimization
initialGeometry = createDMVector(adjointDM);
CHKERR mField.getInterface<VecManager>()->setOtherLocalGhostVector(
"ADJOINT", "ADJOINT_FIELD", "GEOMETRY", RowColData::ROW, initialGeometry,
INSERT_VALUES, SCATTER_FORWARD);
// Setup topology optimization modes
/**
* Generate modes for topology optimization design parameterization
*
* These modes represent perturbations to the geometry that can be used
* as design variables in topology optimization. The modes are defined
* through the Python interface and provide spatial basis functions
* for shape modifications.
*/
CHKERR topologyModes();
// Initialize optimization variables and run solver
/**
* Start optimization with zero initial guess for design variables
*
* The TAO solver will iteratively:
* 1. Evaluate objective function at current design point
* 2. Compute gradients using adjoint sensitivity analysis
* 3. Update design variables using L-BFGS algorithm
* 4. Check convergence criteria
*
* Each function evaluation involves solving the forward elastic problem
* and the adjoint problem to compute sensitivities efficiently.
*/
auto x0 = vectorDuplicate(g);
CHKERR VecSet(x0, 0.0);
CHKERR TaoSetSolution(tao, x0);
CHKERR TaoSetFromOptions(tao);
auto testing = mField.getInterface<TestingInterface>();
CHKERR testing->addTestCommandLine(
"", 1, "gradient_fd",
[this, tao, x0, g, rank, eval_objective_and_gradient]() mutable {
MoFEMFunctionBegin;
PetscReal eps = 1e-4, design = 0;
CHKERR PetscOptionsGetReal(nullptr, nullptr, "-gradient_fd_eps", &eps,
nullptr);
CHKERR PetscOptionsGetReal(nullptr, nullptr, "-gradient_fd_design",
&design, nullptr);
CHKERR VecSet(x0, design);
PetscReal value;
CHKERR eval_objective_and_gradient(tao, x0, &value, g, this);
auto test_x =
createGhostVector(mField.get_comm(), (!rank) ? modeVecs.size() : 0,
modeVecs.size(), 0, nullptr);
CHKERR VecCopy(x0, test_x);
auto direction = vectorDuplicate(test_x);
auto work_gradient = vectorDuplicate(g);
auto evaluate = [&](Vec x, PetscReal &f) {
return eval_objective_and_gradient(tao, x, &f, work_gradient, this);
};
auto tester = mField.getInterface<OperatorsTester>();
for (PetscInt mode = 0; mode < modeVecs.size(); ++mode) {
CHKERR VecZeroEntries(direction);
if (!rank)
CHKERR VecSetValue(direction, mode, 1., INSERT_VALUES);
CHKERR VecAssemblyBegin(direction);
CHKERR VecAssemblyEnd(direction);
PetscReal fd, analytical;
CHKERR tester->checkScalarCentralFiniteDifference(
test_x, direction, g, eps, evaluate, fd, analytical);
MOFEM_LOG("WORLD", Sev::inform)
<< "Shape gradient mode " << mode
<< ": FD=" << std::setprecision(16) << fd
<< " analytical=" << analytical
<< " error=" << std::abs(fd - analytical);
CHKERR mField.getInterface<TestingInterface>()
->checkValue<TestType::DERIVATIVE>("shape gradient mode " +
std::to_string(mode),
analytical, fd, {1e-6, 1e-6});
}
CHKERR eval_objective_and_gradient(tao, x0, &value, g, this);
MoFEMFunctionReturn(0);
});
CHKERR testing->runTests();
if (testing->checkTest("gradient_fd"))
MoFEMFunctionReturnHot(0);
CHKERR TaoSolve(tao);
CHKERR mField.getInterface<TestingInterface>()->checkSolver(
tao, "shape optimisation adjoint");
// Optimization complete - results available in solution vectors
MOFEM_LOG("WORLD", Sev::inform) << "Topology optimization completed";
MoFEMFunctionReturn(0);
}
//! [Run problem]
//! [Read mesh]
/**
* @brief Read mesh from file and setup material/boundary condition meshsets
*
* This function loads the finite element mesh from file and processes
* associated meshsets that define material properties and boundary conditions.
* The mesh is typically generated using CUBIT and exported in .h5m format.
*
* Meshsets are used to group elements/faces by:
* - Material properties (for different material blocks)
* - Boundary conditions (for applying loads and constraints)
* - Optimization regions (for topology optimization)
*
* @return MoFEMErrorCode Success or error code
*/
MoFEMErrorCode Example::readMesh() {
MoFEMFunctionBegin;
auto simple = mField.getInterface<Simple>();
CHKERR simple->getOptions(); // Read command line options
CHKERR simple->loadFile(); // Load mesh from file
// Add meshsets if config file provided
CHKERR mField.getInterface<MeshsetsManager>()->setMeshsetFromFile();
MoFEMFunctionReturn(0);
}
//! [Read mesh]
//! [Set up problem]
/**
* @brief Setup finite element fields, approximation spaces and degrees of freedom
*
* This function configures the finite element problem by:
* 1. Setting up the displacement field "U" with vector approximation
* 2. Setting up the geometry field "GEOMETRY" for mesh deformation
* 3. Defining polynomial approximation order and basis functions
* 4. Creating degrees of freedom on mesh entities
*
* The displacement field uses H1 vector space for standard elasticity.
* The geometry field allows mesh modification during topology optimization.
* Different basis functions (Ainsworth-Legendre vs Demkowicz) can be selected.
*
* @return MoFEMErrorCode Success or error code
*/
MoFEMErrorCode Example::setupProblem() {
MoFEMFunctionBegin;
Simple *simple = mField.getInterface<Simple>();
// Select basis functions for finite element approximation
enum bases { AINSWORTH, DEMKOWICZ, LASBASETOPT };
const char *list_bases[LASBASETOPT] = {"ainsworth", "demkowicz"};
PetscInt choice_base_value = AINSWORTH;
CHKERR PetscOptionsGetEList(PETSC_NULLPTR, NULL, "-base", list_bases,
LASBASETOPT, &choice_base_value, PETSC_NULLPTR);
switch (choice_base_value) {
case AINSWORTH:
base = AINSWORTH_LEGENDRE_BASE;
MOFEM_LOG("WORLD", Sev::inform)
<< "Set AINSWORTH_LEGENDRE_BASE for displacements";
break;
case DEMKOWICZ:
base = DEMKOWICZ_JACOBI_BASE;
MOFEM_LOG("WORLD", Sev::inform)
<< "Set DEMKOWICZ_JACOBI_BASE for displacements";
break;
default:
base = LASTBASE;
break;
}
// Add finite element fields
/**
* Setup displacement field "U" - the primary unknown in elasticity
* This field represents displacement vector at each node/DOF
*/
CHKERR simple->addDomainField("U", H1, base, SPACE_DIM);
CHKERR simple->addBoundaryField("U", H1, base, SPACE_DIM);
CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-order", &fieldOrder,
PETSC_NULLPTR);
/**
* Setup geometry field "GEOMETRY" - used for mesh deformation in optimization
* This field stores current nodal coordinates and can be modified
* during topology optimization to represent design changes
*/
CHKERR simple->addDataField("GEOMETRY", H1, base, SPACE_DIM);
// Set polynomial approximation order for both fields
CHKERR simple->setFieldOrder("U", fieldOrder);
CHKERR simple->setFieldOrder("GEOMETRY", fieldOrder);
CHKERR simple->setUp();
// Project higher-order geometry representation onto geometry field
/**
* For higher-order elements, this projects the exact geometry
* onto the geometry field to maintain curved boundaries accurately
*/
auto project_ho_geometry = [&]() {
Projection10NodeCoordsOnField ent_method(mField, "GEOMETRY");
return mField.loop_dofs("GEOMETRY", ent_method);
};
CHKERR project_ho_geometry();
// Check if plane strain assumption should be used
CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-plane_strain",
&is_plane_strain, PETSC_NULLPTR);
MoFEMFunctionReturn(0);
}
//! [Set up problem]
//! [Setup adjoint]
/**
* @brief Setup adjoint fields and finite elements for sensitivity analysis
*
* The adjoint method is used to efficiently compute gradients of the objective
* function with respect to design variables. This function sets up:
*
* 1. ADJOINT_FIELD - stores adjoint variables (Lagrange multipliers)
* 2. ADJOINT_DM - data manager for adjoint problem
* 3. Adjoint finite elements for domain and boundary
*
* The adjoint equation is: K^T * λ = ∂f/∂u
* where λ are adjoint variables, K is stiffness matrix, f is objective
*
* @return MoFEMErrorCode Success or error code
*/
MoFEMErrorCode Example::setupAdJoint() {
MoFEMFunctionBegin;
auto simple = mField.getInterface<Simple>();
// Create adjoint data manager and field
auto create_adjoint_dm = [&]() {
auto adjoint_dm = createDM(mField.get_comm(), "DMMOFEM");
auto add_field = [&]() {
MoFEMFunctionBegin;
CHKERR mField.add_field("ADJOINT_FIELD", H1, base, SPACE_DIM);
CHKERR mField.add_ents_to_field_by_dim(simple->getMeshset(), SPACE_DIM,
"ADJOINT_FIELD");
for (auto tt = MBEDGE; tt <= moab::CN::TypeDimensionMap[SPACE_DIM].second;
++tt)
CHKERR mField.set_field_order(simple->getMeshset(), tt, "ADJOINT_FIELD",
fieldOrder);
CHKERR mField.set_field_order(simple->getMeshset(), MBVERTEX,
"ADJOINT_FIELD", 1);
CHKERR mField.build_fields();
MoFEMFunctionReturn(0);
};
auto add_adjoint_fe_impl = [&]() {
MoFEMFunctionBegin;
CHKERR mField.add_finite_element("ADJOINT_DOMAIN_FE");
CHKERR mField.modify_finite_element_add_field_row("ADJOINT_DOMAIN_FE",
"ADJOINT_FIELD");
CHKERR mField.modify_finite_element_add_field_col("ADJOINT_DOMAIN_FE",
"ADJOINT_FIELD");
CHKERR mField.modify_finite_element_add_field_data("ADJOINT_DOMAIN_FE",
"ADJOINT_FIELD");
CHKERR mField.modify_finite_element_add_field_data("ADJOINT_DOMAIN_FE",
"GEOMETRY");
CHKERR mField.add_ents_to_finite_element_by_dim(
simple->getMeshset(), SPACE_DIM, "ADJOINT_DOMAIN_FE");
CHKERR mField.build_finite_elements("ADJOINT_DOMAIN_FE");
CHKERR mField.add_finite_element("ADJOINT_BOUNDARY_FE");
CHKERR mField.modify_finite_element_add_field_row("ADJOINT_BOUNDARY_FE",
"ADJOINT_FIELD");
CHKERR mField.modify_finite_element_add_field_col("ADJOINT_BOUNDARY_FE",
"ADJOINT_FIELD");
CHKERR mField.modify_finite_element_add_field_data("ADJOINT_BOUNDARY_FE",
"ADJOINT_FIELD");
CHKERR mField.modify_finite_element_add_field_data("ADJOINT_BOUNDARY_FE",
"GEOMETRY");
auto block_name = "OPTIMISE";
auto mesh_mng = mField.getInterface<MeshsetsManager>();
auto bcs = mesh_mng->getCubitMeshsetPtr(
std::regex((boost::format("%s(.*)") % block_name).str())
);
for (auto bc : bcs) {
CHKERR mField.add_ents_to_finite_element_by_dim(
bc->getMeshset(), SPACE_DIM - 1, "ADJOINT_BOUNDARY_FE");
}
CHKERR mField.build_finite_elements("ADJOINT_BOUNDARY_FE");
CHKERR mField.build_adjacencies(simple->getBitRefLevel(),
simple->getBitRefLevelMask());
MoFEMFunctionReturn(0);
};
auto set_adjoint_dm_imp = [&]() {
MoFEMFunctionBegin;
CHKERR DMMoFEMCreateMoFEM(adjoint_dm, &mField, "ADJOINT",
simple->getBitRefLevel(),
simple->getBitRefLevelMask());
CHKERR DMMoFEMSetDestroyProblem(adjoint_dm, PETSC_TRUE);
CHKERR DMSetFromOptions(adjoint_dm);
CHKERR DMMoFEMAddElement(adjoint_dm, "ADJOINT_DOMAIN_FE");
CHKERR DMMoFEMAddElement(adjoint_dm, "ADJOINT_BOUNDARY_FE");
CHKERR DMMoFEMSetSquareProblem(adjoint_dm, PETSC_TRUE);
CHKERR DMMoFEMSetIsPartitioned(adjoint_dm, PETSC_TRUE);
mField.getInterface<ProblemsManager>()->buildProblemFromFields =
PETSC_TRUE;
CHKERR DMSetUp(adjoint_dm);
mField.getInterface<ProblemsManager>()->buildProblemFromFields =
PETSC_FALSE;
MoFEMFunctionReturn(0);
};
CHK_THROW_MESSAGE(add_field(), "add adjoint field");
CHK_THROW_MESSAGE(add_adjoint_fe_impl(), "add adjoint fe");
CHK_THROW_MESSAGE(set_adjoint_dm_imp(), "set adjoint dm");
return adjoint_dm;
};
adjointDM = create_adjoint_dm();
MoFEMFunctionReturn(0);
}
//
//! [Boundary condition]
MoFEMErrorCode Example::boundaryCondition() {
MoFEMFunctionBegin;
auto simple = mField.getInterface<Simple>();
auto bc_mng = mField.getInterface<BcManager>();
CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_X",
"U", 0, 0);
CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_Y",
"U", 1, 1);
CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_Z",
"U", 2, 2);
CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(),
"REMOVE_ALL", "U", 0, 3);
CHKERR bc_mng->pushMarkDOFsOnEntities<DisplacementCubitBcData>(
simple->getProblemName(), "U");
MoFEMFunctionReturn(0);
}
//! [Boundary condition]
//! [Adjoint modes]
MoFEMErrorCode Example::topologyModes() {
MoFEMFunctionBegin;
auto opt_ents = get_range_from_block(mField, "OPTIMISE", SPACE_DIM - 1);
auto subset_dm_bdy = createDM(mField.get_comm(), "DMMOFEM");
CHKERR DMMoFEMSetSquareProblem(subset_dm_bdy, PETSC_TRUE);
CHKERR DMMoFEMCreateSubDM(subset_dm_bdy, adjointDM, "SUBSET_BDY");
CHKERR DMMoFEMAddElement(subset_dm_bdy, "ADJOINT_BOUNDARY_FE");
CHKERR DMMoFEMAddSubFieldRow(subset_dm_bdy, "ADJOINT_FIELD",
boost::make_shared<Range>(opt_ents));
CHKERR DMMoFEMAddSubFieldCol(subset_dm_bdy, "ADJOINT_FIELD",
boost::make_shared<Range>(opt_ents));
CHKERR DMSetUp(subset_dm_bdy);
auto subset_dm_domain = createDM(mField.get_comm(), "DMMOFEM");
CHKERR DMMoFEMSetSquareProblem(subset_dm_domain, PETSC_TRUE);
CHKERR DMMoFEMCreateSubDM(subset_dm_domain, adjointDM, "SUBSET_DOMAIN");
CHKERR DMMoFEMAddElement(subset_dm_domain, "ADJOINT_DOMAIN_FE");
CHKERR DMMoFEMAddSubFieldRow(subset_dm_domain, "ADJOINT_FIELD");
CHKERR DMMoFEMAddSubFieldCol(subset_dm_domain, "ADJOINT_FIELD");
CHKERR DMSetUp(subset_dm_domain);
// remove dofs on boundary of the domain
auto remove_dofs = [&]() {
MoFEMFunctionBegin;
std::array<Range, 3> remove_dim_ents;
remove_dim_ents[0] =
get_range_from_block(mField, "OPT_REMOVE_X", SPACE_DIM - 1);
remove_dim_ents[1] =
get_range_from_block(mField, "OPT_REMOVE_Y", SPACE_DIM - 1);
remove_dim_ents[2] =
get_range_from_block(mField, "OPT_REMOVE_Z", SPACE_DIM - 1);
for (int d = 0; d != 3; ++d) {
MOFEM_LOG("WORLD", Sev::inform)
<< "Removing topology modes on block OPT_REMOVE_" << (char)('X' + d)
<< " with " << remove_dim_ents[d].size() << " entities";
}
Range body_ents;
CHKERR mField.get_moab().get_entities_by_dimension(0, SPACE_DIM, body_ents,
true);
auto skin = moab::Skinner(&mField.get_moab());
Range boundary_ents;
CHKERR skin.find_skin(0, body_ents, false, boundary_ents);
for (int d = 0; d != 3; ++d) {
boundary_ents = subtract(boundary_ents, remove_dim_ents[d]);
}
ParallelComm *pcomm =
ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
CHKERR pcomm->filter_pstatus(boundary_ents,
PSTATUS_SHARED | PSTATUS_MULTISHARED,
PSTATUS_NOT, -1, nullptr);
for (auto d = SPACE_DIM - 2; d >= 0; --d) {
if (d >= 0) {
Range ents;
CHKERR mField.get_moab().get_adjacencies(boundary_ents, d, false, ents,
moab::Interface::UNION);
boundary_ents.merge(ents);
} else {
Range verts;
CHKERR mField.get_moab().get_connectivity(boundary_ents, verts);
boundary_ents.merge(verts);
}
CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
boundary_ents);
}
boundary_ents.merge(opt_ents);
CHKERR mField.getInterface<ProblemsManager>()->removeDofsOnEntities(
"SUBSET_DOMAIN", "ADJOINT_FIELD", boundary_ents);
for (int d = 0; d != 3; ++d) {
CHKERR mField.getInterface<ProblemsManager>()->removeDofsOnEntities(
"SUBSET_DOMAIN", "ADJOINT_FIELD", remove_dim_ents[d], d, d);
}
// #ifndef NDEBUG
if (mField.get_comm_rank() == 0) {
CHKERR save_range(mField.get_moab(), "topoMode_boundary_ents.vtk",
boundary_ents);
}
// #endif
MoFEMFunctionReturn(0);
};
CHKERR remove_dofs();
auto get_lhs_fe = [&]() {
auto fe_lhs = boost::make_shared<BoundaryEle>(mField);
fe_lhs->getRuleHook = [](int, int, int p_data) {
return 2 * p_data + p_data - 1;
};
auto &pip = fe_lhs->getOpPtrVector();
CHKERR AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(pip, {},
"GEOMETRY");
using OpMass = FormsIntegrators<BoundaryEleOp>::Assembly<A>::BiLinearForm<
I>::OpMass<1, SPACE_DIM>;
pip.push_back(new OpMass("ADJOINT_FIELD", "ADJOINT_FIELD",
[](double, double, double) { return 1.; }));
return fe_lhs;
};
auto get_rhs_fe = [&]() {
auto fe_rhs = boost::make_shared<BoundaryEle>(mField);
fe_rhs->getRuleHook = [](int, int, int p_data) {
return 2 * p_data + p_data - 1;
};
auto &pip = fe_rhs->getOpPtrVector();
CHKERR AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(pip, {},
"GEOMETRY");
return fe_rhs;
};
auto block_name = "OPTIMISE";
auto mesh_mng = mField.getInterface<MeshsetsManager>();
auto bcs = mesh_mng->getCubitMeshsetPtr(
std::regex((boost::format("%s(.*)") % block_name).str())
);
for (auto &v : modeVecs) {
v = createDMVector(subset_dm_bdy);
}
using OP = FormsIntegrators<BoundaryEleOp>::Assembly<A>::OpBase;
struct OpMode : public OP {
OpMode(const std::string name,
boost::shared_ptr<ObjectiveFunctionData> python_ptr, int id,
std::vector<SmartPetscObj<Vec>> mode_vecs,
std::vector<std::array<double, 3>> mode_centroids,
std::vector<std::array<double, 6>> mode_bboxes, int block_counter,
int mode_counter, boost::shared_ptr<Range> range = nullptr)
: OP(name, name, OP::OPROW, range), pythonPtr(python_ptr), iD(id),
modeVecs(mode_vecs), modeCentroids(mode_centroids),
modeBboxes(mode_bboxes), blockCounter(block_counter),
modeCounter(mode_counter) {}
MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
MoFEMFunctionBegin;
if (OP::entsPtr) {
if (OP::entsPtr->find(this->getFEEntityHandle()) == OP::entsPtr->end())
MoFEMFunctionReturnHot(0);
}
auto nb_rows = data.getIndices().size();
if (!nb_rows) {
MoFEMFunctionReturnHot(0);
}
auto nb_base_functions = data.getN().size2();
FTENSOR_INDEX(SPACE_DIM, i);
CHKERR pythonPtr->blockModes(iD, OP::getCoordsAtGaussPts(),
modeCentroids[blockCounter],
modeBboxes[blockCounter], blockModes);
auto nb_integration_pts = getGaussPts().size2();
if (blockModes.size2() != 3 * nb_integration_pts) {
MOFEM_LOG("WORLD", Sev::error)
<< "Number of modes does not match number of integration points: "
<< blockModes.size2() << "!=" << 3 * nb_integration_pts;
CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY, "modes/integration points");
}
VectorDouble nf(nb_rows);
int nb_modes = blockModes.size1();
for (auto mode = 0; mode != nb_modes; ++mode) {
nf.clear();
// get mode
auto t_mode = getFTensor1FromPtr<3>(&blockModes(mode, 0));
// get element volume
const double vol = OP::getMeasure();
// get integration weights
auto t_w = OP::getFTensor0IntegrationWeight();
// get base function gradient on rows
auto t_base = data.getFTensor0N();
// loop over integration points
for (int gg = 0; gg != nb_integration_pts; gg++) {
// take into account Jacobian
const double alpha = t_w * vol;
// loop over rows base functions
auto t_nf = getFTensor1FromPtr<SPACE_DIM>(nf.data().data());
int rr = 0;
for (; rr != nb_rows / SPACE_DIM; ++rr) {
t_nf(i) += alpha * t_base * t_mode(i);
++t_base;
++t_nf;
}
for (; rr < nb_base_functions; ++rr)
++t_base;
++t_w; // move to another integration weight
++t_mode; // move to another mode
}
Vec vec = modeVecs[modeCounter + mode];
auto size = data.getIndices().size();
auto *indices = data.getIndices().data().data();
auto *nf_data = nf.data().data();
CHKERR VecSetValues(vec, size, indices, nf_data, ADD_VALUES);
}
MoFEMFunctionReturn(0);
}
private:
boost::shared_ptr<ObjectiveFunctionData> pythonPtr;
MatrixDouble blockModes;
std::vector<std::array<double, 3>> modeCentroids;
std::vector<std::array<double, 6>> modeBboxes;
int iD;
std::vector<SmartPetscObj<Vec>> modeVecs;
int blockCounter;
int modeCounter;
};
auto solve_bdy = [&]() {
MoFEMFunctionBegin;
auto fe_lhs = get_lhs_fe();
auto fe_rhs = get_rhs_fe();
int block_counter = 0;
int mode_counter = 0;
for (auto &bc : bcs) {
auto id = bc->getMeshsetId();
Range ents;
CHKERR mField.get_moab().get_entities_by_handle(bc->getMeshset(), ents,
true);
auto range = boost::make_shared<Range>(ents);
auto &pip_rhs = fe_rhs->getOpPtrVector();
pip_rhs.push_back(new OpMode("ADJOINT_FIELD", pythonPtr, id, modeVecs,
modeCentroids, modeBBoxes, block_counter,
mode_counter, range));
CHKERR DMoFEMLoopFiniteElements(subset_dm_bdy, "ADJOINT_BOUNDARY_FE",
fe_rhs);
pip_rhs.pop_back();
int nb_modes;
CHKERR pythonPtr->numberOfModes(id, nb_modes);
++block_counter;
mode_counter += nb_modes;
MOFEM_LOG("WORLD", Sev::inform)
<< "Setting mode block block: " << bc->getName()
<< " with ID: " << bc->getMeshsetId()
<< " total modes: " << mode_counter;
}
for (auto &v : modeVecs) {
CHKERR VecAssemblyBegin(v);
CHKERR VecAssemblyEnd(v);
CHKERR VecGhostUpdateBegin(v, ADD_VALUES, SCATTER_REVERSE);
CHKERR VecGhostUpdateEnd(v, ADD_VALUES, SCATTER_REVERSE);
}
auto M = createDMMatrix(subset_dm_bdy);
fe_lhs->B = M;
CHKERR DMoFEMLoopFiniteElements(subset_dm_bdy, "ADJOINT_BOUNDARY_FE",
fe_lhs);
CHKERR MatAssemblyBegin(M, MAT_FINAL_ASSEMBLY);
CHKERR MatAssemblyEnd(M, MAT_FINAL_ASSEMBLY);
auto solver = createKSP(mField.get_comm());
CHKERR KSPSetOperators(solver, M, M);
CHKERR KSPSetFromOptions(solver);
CHKERR KSPSetUp(solver);
auto v = createDMVector(subset_dm_bdy);
for (auto &f : modeVecs) {
CHKERR KSPSolve(solver, f, v);
CHKERR VecSwap(f, v);
}
for (auto &v : modeVecs) {
CHKERR VecGhostUpdateBegin(v, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(v, INSERT_VALUES, SCATTER_FORWARD);
}
MoFEMFunctionReturn(0);
};
CHKERR solve_bdy();
auto get_elastic_fe_lhs = [&]() {
auto fe = boost::make_shared<DomainEle>(mField);
fe->getRuleHook = [](int, int, int p_data) {
return 2 * p_data + p_data - 1;
};
auto &pip = fe->getOpPtrVector();
CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {H1},
"GEOMETRY");
CHKERR HookeOps::opFactoryDomainLhs<SPACE_DIM, A, I, DomainEleOp>(
mField, pip, "ADJOINT_FIELD", "MAT_ADJOINT", Sev::noisy);
return fe;
};
auto get_elastic_fe_rhs = [&]() {
auto fe = boost::make_shared<DomainEle>(mField);
fe->getRuleHook = [](int, int, int p_data) {
return 2 * p_data + p_data - 1;
};
auto &pip = fe->getOpPtrVector();
CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {H1},
"GEOMETRY");
CHKERR HookeOps::opFactoryDomainRhs<SPACE_DIM, A, I, DomainEleOp>(
mField, pip, "ADJOINT_FIELD", "MAT_ADJOINT", Sev::noisy);
return fe;
};
auto adjoint_gradient_postprocess = [&](auto mode) {
MoFEMFunctionBegin;
auto post_proc_mesh = boost::make_shared<moab::Core>();
auto post_proc_begin =
boost::make_shared<PostProcBrokenMeshInMoabBaseBegin>(mField,
post_proc_mesh);
auto post_proc_end = boost::make_shared<PostProcBrokenMeshInMoabBaseEnd>(
mField, post_proc_mesh);
auto geom_vec = boost::make_shared<MatrixDouble>();
auto post_proc_fe =
boost::make_shared<PostProcEleDomain>(mField, post_proc_mesh);
AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
post_proc_fe->getOpPtrVector(), {H1}, "GEOMETRY");
post_proc_fe->getOpPtrVector().push_back(
new OpCalculateVectorFieldValues<SPACE_DIM>("ADJOINT_FIELD", geom_vec,
modeVecs[mode]));
using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
post_proc_fe->getOpPtrVector().push_back(
new OpPPMap(
post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
{},
{{"MODE", geom_vec}},
{},
{}
)
);
CHKERR DMoFEMPreProcessFiniteElements(adjointDM,
post_proc_begin->getFEMethod());
CHKERR DMoFEMLoopFiniteElements(adjointDM, "ADJOINT_DOMAIN_FE",
post_proc_fe);
CHKERR DMoFEMPostProcessFiniteElements(adjointDM,
post_proc_begin->getFEMethod());
CHKERR post_proc_end->writeFile("mode_" + std::to_string(mode) + ".h5m");
MoFEMFunctionReturn(0);
};
auto solve_domain = [&]() {
MoFEMFunctionBegin;
auto fe_lhs = get_elastic_fe_lhs();
auto fe_rhs = get_elastic_fe_rhs();
auto v = createDMVector(subset_dm_domain);
auto F = vectorDuplicate(v);
fe_rhs->f = F;
auto M = createDMMatrix(subset_dm_domain);
fe_lhs->B = M;
CHKERR DMoFEMLoopFiniteElements(subset_dm_domain, "ADJOINT_DOMAIN_FE",
fe_lhs);
CHKERR MatAssemblyBegin(M, MAT_FINAL_ASSEMBLY);
CHKERR MatAssemblyEnd(M, MAT_FINAL_ASSEMBLY);
auto solver = createKSP(mField.get_comm());
CHKERR KSPSetOperators(solver, M, M);
CHKERR KSPSetFromOptions(solver);
CHKERR KSPSetUp(solver);
int mode_counter = 0;
for (auto &f : modeVecs) {
CHKERR mField.getInterface<FieldBlas>()->setField(0, "ADJOINT_FIELD");
CHKERR DMoFEMMeshToLocalVector(subset_dm_bdy, f, INSERT_VALUES,
SCATTER_REVERSE);
CHKERR VecZeroEntries(F);
CHKERR DMoFEMLoopFiniteElements(subset_dm_domain, "ADJOINT_DOMAIN_FE",
fe_rhs);
CHKERR VecAssemblyBegin(F);
CHKERR VecAssemblyEnd(F);
CHKERR VecGhostUpdateBegin(F, ADD_VALUES, SCATTER_REVERSE);
CHKERR VecGhostUpdateEnd(F, ADD_VALUES, SCATTER_REVERSE);
CHKERR KSPSolve(solver, F, v);
CHKERR VecGhostUpdateBegin(v, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(v, INSERT_VALUES, SCATTER_FORWARD);
CHKERR DMoFEMMeshToLocalVector(subset_dm_domain, v, INSERT_VALUES,
SCATTER_REVERSE);
auto m = createDMVector(adjointDM);
CHKERR DMoFEMMeshToLocalVector(adjointDM, m, INSERT_VALUES,
SCATTER_FORWARD);
f = m;
++mode_counter;
}
MoFEMFunctionReturn(0);
};
CHKERR solve_domain();
for (int i = 0; i < modeVecs.size(); ++i) {
CHKERR adjoint_gradient_postprocess(i);
}
MoFEMFunctionReturn(0);
}
//! [Adjoint modes]
//! [Push operators to pipeline]
MoFEMErrorCode Example::assembleSystem() {
MoFEMFunctionBegin;
auto pip = mField.getInterface<PipelineManager>();
//! [Integration rule]
CHKERR pip->setDomainRhsIntegrationRule(integrationRule);
CHKERR pip->setDomainLhsIntegrationRule(integrationRule);
CHKERR pip->setBoundaryRhsIntegrationRule(integrationRule);
CHKERR pip->setBoundaryLhsIntegrationRule(integrationRule);
//! [Integration rule]
CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
pip->getOpDomainLhsPipeline(), {H1}, "GEOMETRY");
CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
pip->getOpDomainRhsPipeline(), {H1}, "GEOMETRY");
CHKERR AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
pip->getOpBoundaryRhsPipeline(), {NOSPACE}, "GEOMETRY");
CHKERR AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
pip->getOpBoundaryLhsPipeline(), {NOSPACE}, "GEOMETRY");
//! [Push domain stiffness matrix to pipeline]
// Add LHS operator for elasticity (stiffness matrix)
CHKERR HookeOps::opFactoryDomainLhs<SPACE_DIM, A, I, DomainEleOp>(
mField, pip->getOpDomainLhsPipeline(), "U", "MAT_ELASTIC", Sev::noisy);
//! [Push domain stiffness matrix to pipeline]
// Add RHS operator for internal forces
CHKERR HookeOps::opFactoryDomainRhs<SPACE_DIM, A, I, DomainEleOp>(
mField, pip->getOpDomainRhsPipeline(), "U", "MAT_ELASTIC", Sev::noisy);
//! [Push Body forces]
CHKERR DomainRhsBCs::AddFluxToPipeline<OpDomainRhsBCs>::add(
pip->getOpDomainRhsPipeline(), mField, "U", Sev::inform);
//! [Push Body forces]
//! [Push natural boundary conditions]
// Add force boundary condition
CHKERR BoundaryRhsBCs::AddFluxToPipeline<OpBoundaryRhsBCs>::add(
pip->getOpBoundaryRhsPipeline(), mField, "U", -1, Sev::inform);
// Add case for mix type of BCs
CHKERR BoundaryLhsBCs::AddFluxToPipeline<OpBoundaryLhsBCs>::add(
pip->getOpBoundaryLhsPipeline(), mField, "U", Sev::noisy);
//! [Push natural boundary conditions]
MoFEMFunctionReturn(0);
}
//! [Push operators to pipeline]
//! [Solve]
MoFEMErrorCode Example::solveElastic() {
MoFEMFunctionBegin;
auto simple = mField.getInterface<Simple>();
auto dm = simple->getDM();
auto f = createDMVector(dm);
auto d = vectorDuplicate(f);
CHKERR VecZeroEntries(d);
CHKERR DMoFEMMeshToLocalVector(dm, d, INSERT_VALUES, SCATTER_REVERSE);
auto set_essential_bc = [&]() {
MoFEMFunctionBegin;
// This is low level pushing finite elements (pipelines) to solver
auto ksp_ctx_ptr = getDMKspCtx(dm);
auto pre_proc_rhs = boost::make_shared<FEMethod>();
auto post_proc_rhs = boost::make_shared<FEMethod>();
auto post_proc_lhs = boost::make_shared<FEMethod>();
auto get_pre_proc_hook = [&]() {
return EssentialPreProc<DisplacementCubitBcData>(mField, pre_proc_rhs,
{});
};
pre_proc_rhs->preProcessHook = get_pre_proc_hook();
auto get_post_proc_hook_rhs = [this, post_proc_rhs]() {
MoFEMFunctionBeginHot;
CHKERR EssentialPostProcRhs<DisplacementCubitBcData>(mField,
post_proc_rhs, 1.)();
MoFEMFunctionReturnHot(0);
};
auto get_post_proc_hook_lhs = [this, post_proc_lhs]() {
MoFEMFunctionBeginHot;
CHKERR EssentialPostProcLhs<DisplacementCubitBcData>(mField,
post_proc_lhs, 1.)();
MoFEMFunctionReturnHot(0);
};
post_proc_rhs->postProcessHook = get_post_proc_hook_rhs;
post_proc_lhs->postProcessHook = get_post_proc_hook_lhs;
ksp_ctx_ptr->getPreProcComputeRhs().push_front(pre_proc_rhs);
ksp_ctx_ptr->getPostProcComputeRhs().push_back(post_proc_rhs);
ksp_ctx_ptr->getPostProcSetOperators().push_back(post_proc_lhs);
MoFEMFunctionReturn(0);
};
auto setup_and_solve = [&](auto solver) {
MoFEMFunctionBegin;
BOOST_LOG_SCOPED_THREAD_ATTR("Timeline", attrs::timer());
MOFEM_LOG("TIMER", Sev::noisy) << "KSPSetUp";
CHKERR KSPSetUp(solver);
MOFEM_LOG("TIMER", Sev::noisy) << "KSPSetUp <= Done";
MOFEM_LOG("TIMER", Sev::noisy) << "KSPSolve";
CHKERR KSPSolve(solver, f, d);
MOFEM_LOG("TIMER", Sev::noisy) << "KSPSolve <= Done";
MoFEMFunctionReturn(0);
};
MOFEM_LOG_CHANNEL("TIMER");
MOFEM_LOG_TAG("TIMER", "timer");
CHKERR set_essential_bc();
CHKERR setup_and_solve(kspElastic);
CHKERR VecGhostUpdateBegin(d, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(d, INSERT_VALUES, SCATTER_FORWARD);
CHKERR DMoFEMMeshToLocalVector(dm, d, INSERT_VALUES, SCATTER_REVERSE);
auto evaluate_field_at_the_point = [&]() {
MoFEMFunctionBegin;
int coords_dim = 3;
std::array<double, 3> field_eval_coords{0.0, 0.0, 0.0};
PetscBool do_eval_field = PETSC_FALSE;
CHKERR PetscOptionsGetRealArray(NULL, NULL, "-field_eval_coords",
field_eval_coords.data(), &coords_dim,
&do_eval_field);
if (do_eval_field) {
auto vector_field_ptr = boost::make_shared<MatrixDouble>();
auto field_eval_data =
mField.getInterface<FieldEvaluatorInterface>()->getData<DomainEle>();
CHKERR mField.getInterface<FieldEvaluatorInterface>()
->buildTree<SPACE_DIM>(field_eval_data, simple->getDomainFEName());
field_eval_data->setEvalPoints(field_eval_coords.data(), 1);
auto no_rule = [](int, int, int) { return -1; };
auto field_eval_fe_ptr = field_eval_data->feMethodPtr;
field_eval_fe_ptr->getRuleHook = no_rule;
field_eval_fe_ptr->getOpPtrVector().push_back(
new OpCalculateVectorFieldValues<SPACE_DIM>("U", vector_field_ptr));
CHKERR mField.getInterface<FieldEvaluatorInterface>()
->evalFEAtThePoint<SPACE_DIM>(
field_eval_coords.data(), 1e-12, simple->getProblemName(),
simple->getDomainFEName(), field_eval_data,
mField.get_comm_rank(), mField.get_comm_rank(), nullptr, MF_EXIST,
QUIET);
if (vector_field_ptr->size1()) {
auto t_disp = getFTensor1FromMat<SPACE_DIM>(*vector_field_ptr);
if constexpr (SPACE_DIM == 2)
MOFEM_LOG("FieldEvaluator", Sev::inform)
<< "U_X: " << t_disp(0) << " U_Y: " << t_disp(1);
else
MOFEM_LOG("FieldEvaluator", Sev::inform)
<< "U_X: " << t_disp(0) << " U_Y: " << t_disp(1)
<< " U_Z: " << t_disp(2);
}
MOFEM_LOG_SYNCHRONISE(mField.get_comm());
}
MoFEMFunctionReturn(0);
};
CHKERR evaluate_field_at_the_point();
MoFEMFunctionReturn(0);
}
//! [Solve]
//! [Postprocess results]
MoFEMErrorCode Example::postprocessElastic(int iter,
SmartPetscObj<Vec> adjoint_vector) {
MoFEMFunctionBegin;
auto simple = mField.getInterface<Simple>();
CHKERR ElasticPostProc::postProcessElasticResults<SPACE_DIM, DomainEle,
BoundaryEle, SideEle>(
mField, simple->getDM(), simple->getDomainFEName(),
"out_elastic_" + std::to_string(iter) + ".h5m",
{{"ADJOINT", adjoint_vector}}, {}, Sev::noisy);
MoFEMFunctionReturn(0);
}
//! [Postprocess results]
//! [calculateGradient]
using DomainBaseOp = FormsIntegrators<DomainEleOp>::Assembly<A>::OpBase;
template <int SPACE_DIM>
struct OpAdJointGradTimesSymTensor<SPACE_DIM, IntegrationType::GAUSS,
DomainBaseOp> : public DomainBaseOp {
using OP = DomainBaseOp;
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)
: OP(field_name, field_name, DomainEleOp::OPROW), commPtr(comm_ptr),
jac(jac), diffJac(diff_jac), cofVals(cof_vals) {}
protected:
boost::shared_ptr<HookeOps::CommonData> commPtr;
boost::shared_ptr<MatrixDouble> jac;
boost::shared_ptr<MatrixDouble> diffJac;
boost::shared_ptr<VectorDouble> cofVals;
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data);
};
template <int SPACE_DIM>
MoFEMErrorCode
OpAdJointGradTimesSymTensor<SPACE_DIM, IntegrationType::GAUSS, DomainBaseOp>::
iNtegrate(EntitiesFieldData::EntData &row_data) {
MoFEMFunctionBegin;
FTENSOR_INDEX(SPACE_DIM, i);
FTENSOR_INDEX(SPACE_DIM, j);
FTENSOR_INDEX(SPACE_DIM, k);
FTENSOR_INDEX(SPACE_DIM, l);
FTENSOR_INDEX(SPACE_DIM, I);
FTENSOR_INDEX(SPACE_DIM, J);
auto t_diff_symm = FTensor::diff_symmetrize<double>();
// get element volume
const double vol = OP::getMeasure();
// get integration weights
auto t_w = OP::getFTensor0IntegrationWeight();
// get Jacobian values
auto t_jac = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*(jac));
// get diff Jacobian values
auto t_diff_jac = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*(diffJac));
// get cofactor values
auto t_cof = getFTensor0FromVec(*(cofVals));
// get base function gradient on rows
auto t_row_grad = row_data.getFTensor1DiffN<SPACE_DIM>();
// get fradient of the field
auto t_grad_u =
getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*(commPtr->matGradPtr));
// get field gradient values
auto t_cauchy_stress =
getFTensor2SymmetricFromMat<SPACE_DIM>(*(commPtr->getMatCauchyStress()));
// material stiffness tensor
auto t_D =
getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(*(commPtr->matDPtr));
// loop over integration points
for (int gg = 0; gg != OP::nbIntegrationPts; gg++) {
// take into account Jacobian
const double alpha = t_w * vol;
auto t_det = determinantTensor(t_jac);
FTensor::Tensor2<double, SPACE_DIM, SPACE_DIM> t_inv_jac;
CHKERR invertTensor(t_jac, t_det, t_inv_jac);
// Calculate the variation of the gradient due to geometry change
FTensor::Tensor2<double, SPACE_DIM, SPACE_DIM> t_diff_inv_jac;
t_diff_inv_jac(i, j) =
-(t_inv_jac(i, I) * t_diff_jac(I, J)) * t_inv_jac(J, j);
FTensor::Tensor2<double, SPACE_DIM, SPACE_DIM> t_diff_grad;
t_diff_grad(i, j) = t_grad_u(i, k) * t_diff_inv_jac(k, j);
// Calculate the variation of the strain tensor
FTensor::Tensor2<double, SPACE_DIM, SPACE_DIM> t_diff_strain;
t_diff_strain(i, j) = t_diff_symm(i, j, k, l) * t_diff_grad(k, l);
// Calculate the variation of the stress tensor
FTensor::Tensor2<double, SPACE_DIM, SPACE_DIM> t_diff_stress;
t_diff_stress(i, j) = t_D(i, j, k, l) * t_diff_strain(k, l);
// get rhs vector
auto t_nf = OP::template getNf<SPACE_DIM>();
// loop over rows base functions
int rr = 0;
for (; rr != OP::nbRows / SPACE_DIM; rr++) {
FTensor::Tensor1<double, SPACE_DIM> t_diff_row_grad;
t_diff_row_grad(k) = t_row_grad(j) * t_diff_inv_jac(j, k);
// Variation due to change in geometry (diff base)
t_nf(j) += alpha * t_diff_row_grad(i) * t_cauchy_stress(i, j);
// Variation due to change in domain (cofactor)
t_nf(j) += (alpha * t_cof) * t_row_grad(i) * t_cauchy_stress(i, j);
// Variation due to change in stress (diff stress)
t_nf(j) += alpha * t_row_grad(i) * t_diff_stress(i, j);
++t_row_grad;
++t_nf;
}
for (; rr < OP::nbRowBaseFunctions; ++rr) {
++t_row_grad;
}
++t_grad_u;
++t_cauchy_stress;
++t_jac;
++t_diff_jac;
++t_cof;
++t_w; // move to another integration weight
}
MoFEMFunctionReturn(0);
}
//J = 1/2 * sigma : epsilon
// calculate dJ/du = (dJ/dsigma * dsigma/depsilon + dJ/depsilon) depsilon / dgrad_u * dgrad_u / du
struct OpStateSensitivity : public DomainBaseOp {
using OP = DomainBaseOp;
OpStateSensitivity(const std::string field_name,
boost::shared_ptr<ObjectiveFunctionData> python_ptr,
boost::shared_ptr<HookeOps::CommonData> comm_ptr,
boost::shared_ptr<MatrixDouble> u_ptr)
: OP(field_name, field_name, DomainEleOp::OPROW), pythonPtr(python_ptr),
commPtr(comm_ptr), uPtr(u_ptr) {}
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data) {
MoFEMFunctionBegin;
FTENSOR_INDEX(SPACE_DIM, i);
FTENSOR_INDEX(SPACE_DIM, j);
FTENSOR_INDEX(SPACE_DIM, k);
FTENSOR_INDEX(SPACE_DIM, l);
constexpr auto symm_size = (SPACE_DIM * (SPACE_DIM + 1)) / 2;
auto nb_gauss_pts = getGaussPts().size2();
auto objective_dstress =
boost::make_shared<MatrixDouble>(nb_gauss_pts, symm_size);
auto objective_dstrain =
boost::make_shared<MatrixDouble>(nb_gauss_pts, symm_size);
auto objective_du =
boost::make_shared<MatrixDouble>(nb_gauss_pts, SPACE_DIM);
auto evaluate_python = [&]() {
MoFEMFunctionBegin;
auto &coords = OP::getCoordsAtGaussPts();
CHKERR pythonPtr->evalInteriorObjectiveGradientStress(
coords, uPtr, commPtr->getMatCauchyStress(), commPtr->getMatStrain(),
objective_dstress);
CHKERR pythonPtr->evalInteriorObjectiveGradientStrain(
coords, uPtr, commPtr->getMatCauchyStress(), commPtr->getMatStrain(),
objective_dstrain);
CHKERR pythonPtr->evalInteriorObjectiveGradientU(
coords, uPtr, commPtr->getMatCauchyStress(), commPtr->getMatStrain(),
objective_du);
auto vol = OP::getMeasure();
auto t_w = OP::getFTensor0IntegrationWeight();
auto t_D =
getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(*(commPtr->matDPtr));
auto t_row_grad = data.getFTensor1DiffN<SPACE_DIM>();
auto t_row_base = data.getFTensor0N();
auto t_obj_dstress =
getFTensor2SymmetricFromMat<SPACE_DIM>(*objective_dstress);
auto t_obj_dstrain =
getFTensor2SymmetricFromMat<SPACE_DIM>(*objective_dstrain);
auto t_obj_du = getFTensor1FromMat<SPACE_DIM>(*objective_du);
for (int gg = 0; gg != nb_gauss_pts; ++gg) {
const double alpha = t_w * vol;
FTensor::Tensor2<double, SPACE_DIM, SPACE_DIM> t_adjoint_stress;
t_adjoint_stress(i, j) =
t_D(i, j, k, l) * t_obj_dstress(k, l) + t_obj_dstrain(i, j);
auto t_nf = OP::template getNf<SPACE_DIM>();
int rr = 0;
for (; rr != OP::nbRows / SPACE_DIM; rr++) {
t_nf(j) += alpha * t_row_grad(i) * t_adjoint_stress(i, j);
t_nf(j) += alpha * t_row_base * t_obj_du(j);
++t_row_grad;
++t_row_base;
++t_nf;
}
for (; rr < OP::nbRowBaseFunctions; ++rr) {
++t_row_grad;
++t_row_base;
}
++t_obj_dstrain;
++t_obj_dstress;
++t_obj_du;
++t_w;
}
MoFEMFunctionReturn(0);
};
CHKERR evaluate_python();
MoFEMFunctionReturn(0);
}
private:
boost::shared_ptr<ObjectiveFunctionData> pythonPtr;
boost::shared_ptr<HookeOps::CommonData> commPtr;
boost::shared_ptr<MatrixDouble> uPtr;
};
struct OpAdJointObjective : public ForcesAndSourcesCore::UserDataOperator {
using OP = ForcesAndSourcesCore::UserDataOperator;
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> d_u_ptr,
boost::shared_ptr<MatrixDouble> u_ptr,
boost::shared_ptr<double> glob_objective_ptr,
boost::shared_ptr<double> glob_objective_grad_ptr)
: OP(NOSPACE, OP::OPSPACE), pythonPtr(python_ptr), commPtr(comm_ptr),
jacPtr(jac_ptr), diffJacPtr(diff_jac), cofVals(cof_vals),
dGradPtr(d_grad_ptr), dUPtr(d_u_ptr), uPtr(u_ptr),
globObjectivePtr(glob_objective_ptr),
globObjectiveGradPtr(glob_objective_grad_ptr) {}
/**
* @brief Compute objective function contributions at element level
*
* Evaluates Python objective function with current displacement and stress
* state, and accumulates global objective value and gradients.
*/
MoFEMErrorCode doWork(int side, EntityType type,
EntitiesFieldData::EntData &data) {
MoFEMFunctionBegin;
// Define tensor indices for calculations
FTENSOR_INDEX(SPACE_DIM, i);
FTENSOR_INDEX(SPACE_DIM, j);
FTENSOR_INDEX(SPACE_DIM, k);
FTENSOR_INDEX(SPACE_DIM, l);
FTENSOR_INDEX(SPACE_DIM, I);
FTENSOR_INDEX(SPACE_DIM, J);
constexpr auto symm_size = (SPACE_DIM * (SPACE_DIM + 1)) /
2; // size of symmetric tensor in Voigt notation
auto t_diff_symm = FTensor::diff_symmetrize<double>();
// fourth order tensor for symmetrization to convert from gradient to strain
auto nb_gauss_pts =
getGaussPts().size2(); // number of gauss points in the element
auto objective_ptr = boost::make_shared<MatrixDouble>(
1, nb_gauss_pts); // objective function values at gauss points
auto objective_dstress = boost::make_shared<MatrixDouble>(
nb_gauss_pts,
symm_size); // objective function gradient w.r.t. stress at gauss points
auto objective_dstrain = boost::make_shared<MatrixDouble>(
nb_gauss_pts,
symm_size); // objective function gradient w.r.t. strain at gauss points
auto objective_du =
boost::make_shared<MatrixDouble>(nb_gauss_pts, SPACE_DIM);
// We have dJ/dsigma and dJ/depsilon at gauss points, we need to convert it to dJ/du
auto evaluate_python = [&]() {
MoFEMFunctionBegin;
auto &coords = OP::getCoordsAtGaussPts();
CHKERR pythonPtr->evalInteriorObjectiveFunction(
coords, uPtr, commPtr->getMatCauchyStress(), commPtr->getMatStrain(),
objective_ptr);
CHKERR pythonPtr->evalInteriorObjectiveGradientStress(
coords, uPtr, commPtr->getMatCauchyStress(), commPtr->getMatStrain(),
objective_dstress);
CHKERR pythonPtr->evalInteriorObjectiveGradientStrain(
coords, uPtr, commPtr->getMatCauchyStress(), commPtr->getMatStrain(),
objective_dstrain);
CHKERR pythonPtr->evalInteriorObjectiveGradientU(
coords, uPtr, commPtr->getMatCauchyStress(), commPtr->getMatStrain(),
objective_du);
auto t_grad_u =
getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*(commPtr->matGradPtr));
auto t_D =
getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(*(commPtr->matDPtr));
auto t_jac = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*(jacPtr));
auto t_diff_jac = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*(diffJacPtr));
auto t_cof = getFTensor0FromVec(*(cofVals));
auto t_d_grad = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*(dGradPtr));
auto t_obj = getFTensor0FromMat(*objective_ptr);
auto t_obj_dstress =
getFTensor2SymmetricFromMat<SPACE_DIM>(*objective_dstress);
auto t_obj_dstrain =
getFTensor2SymmetricFromMat<SPACE_DIM>(*objective_dstrain);
auto t_obj_du = getFTensor1FromMat<SPACE_DIM>(*objective_du);
auto t_d_u = getFTensor1FromMat<SPACE_DIM>(*dUPtr);
auto vol = OP::getMeasure();
auto t_w = getFTensor0IntegrationWeight();
for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
auto t_det = determinantTensor(t_jac);
FTensor::Tensor2<double, SPACE_DIM, SPACE_DIM> t_inv_jac;
CHKERR invertTensor(t_jac, t_det, t_inv_jac);
FTensor::Tensor2<double, SPACE_DIM, SPACE_DIM> t_diff_inv_jac;
t_diff_inv_jac(i, j) =
-(t_inv_jac(i, I) * t_diff_jac(I, J)) * t_inv_jac(J, j);
FTensor::Tensor2<double, SPACE_DIM, SPACE_DIM> t_diff_grad;
t_diff_grad(i, j) = t_grad_u(i, k) * t_diff_inv_jac(k, j);
FTensor::Tensor2<double, SPACE_DIM, SPACE_DIM> t_d_strain;
t_d_strain(i, j) = t_diff_symm(i, j, k, l) * (
t_d_grad(k, l)
+
t_diff_grad(k, l)
);
auto alpha = t_w * vol;
(*globObjectivePtr) += alpha * t_obj;
(*globObjectiveGradPtr) +=
alpha *
(
t_obj_dstress(i, j) * (t_D(i, j, k, l) * t_d_strain(k, l))
+
t_obj_dstrain(i, j) * t_d_strain(i, j)
+
t_obj_du(i) * t_d_u(i)
+
t_obj * t_cof
);
++t_w;
++t_jac;
++t_diff_jac;
++t_cof;
++t_obj;
++t_obj_dstress;
++t_obj_dstrain;
++t_obj_du;
++t_grad_u;
++t_d_grad;
++t_d_u;
}
MoFEMFunctionReturn(0);
};
CHKERR evaluate_python();
MoFEMFunctionReturn(0);
}
private:
boost::shared_ptr<ObjectiveFunctionData> pythonPtr;
boost::shared_ptr<HookeOps::CommonData> commPtr;
boost::shared_ptr<MatrixDouble> jacPtr;
boost::shared_ptr<MatrixDouble> diffJacPtr;
boost::shared_ptr<VectorDouble> cofVals;
boost::shared_ptr<MatrixDouble> dGradPtr;
boost::shared_ptr<MatrixDouble> dUPtr;
boost::shared_ptr<MatrixDouble> uPtr;
boost::shared_ptr<double> globObjectivePtr;
boost::shared_ptr<double> globObjectiveGradPtr;
};
MoFEMErrorCode Example::calculateGradient(PetscReal *objective_function_value,
Vec objective_function_gradient,
Vec adjoint_vector) {
MOFEM_LOG_CHANNEL("WORLD");
MoFEMFunctionBegin;
auto simple = mField.getInterface<Simple>();
auto ents = get_range_from_block(mField, "OPTIMISE", SPACE_DIM - 1);
auto get_essential_fe = [this]() {
auto post_proc_rhs = boost::make_shared<FEMethod>();
auto get_post_proc_hook_rhs = [this, post_proc_rhs]() {
MoFEMFunctionBeginHot;
CHKERR EssentialPostProcRhs<DisplacementCubitBcData>(mField,
post_proc_rhs, 0)();
MoFEMFunctionReturnHot(0);
};
post_proc_rhs->postProcessHook = get_post_proc_hook_rhs;
return post_proc_rhs;
};
auto get_fd_direvative_fe = [&]() {
auto fe = boost::make_shared<DomainEle>(mField);
fe->getRuleHook = integrationRule;
auto &pip = fe->getOpPtrVector();
AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {H1}, "GEOMETRY");
// Add RHS operators for internal forces
HookeOps::opFactoryDomainRhs<SPACE_DIM, A, I, DomainEleOp>(
mField, pip, "U", "MAT_ELASTIC", Sev::noisy);
return fe;
};
auto calulate_fd_residual = [&](auto eps, auto diff_vec, auto fd_vec) {
MoFEMFunctionBegin;
constexpr bool debug = false;
auto geom_norm = [](MoFEM::Interface &mField) {
MoFEMFunctionBegin;
auto field_blas = mField.getInterface<FieldBlas>();
double nrm2 = 0.0;
auto norm2_field = [&](const double val) {
nrm2 += val * val;
return val;
};
CHKERR field_blas->fieldLambdaOnValues(norm2_field, "GEOMETRY");
MPI_Allreduce(MPI_IN_PLACE, &nrm2, 1, MPI_DOUBLE, MPI_SUM,
mField.get_comm());
MOFEM_LOG("WORLD", Sev::inform) << "Geometry norm: " << sqrt(nrm2);
MoFEMFunctionReturn(0);
};
if constexpr (debug)
CHKERR geom_norm(mField);
auto initial_current_geometry = createDMVector(adjointDM);
CHKERR mField.getInterface<VecManager>()->setOtherLocalGhostVector(
"ADJOINT", "ADJOINT_FIELD", "GEOMETRY", RowColData::ROW,
initial_current_geometry, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecAssemblyBegin(initial_current_geometry);
CHKERR VecAssemblyEnd(initial_current_geometry);
if constexpr (debug)
CHKERR geom_norm(mField);
auto perturb_geometry = [&](auto eps, auto diff_vec) {
MoFEMFunctionBegin;
auto current_geometry = vectorDuplicate(initial_current_geometry);
CHKERR VecCopy(initial_current_geometry, current_geometry);
CHKERR VecAXPY(current_geometry, eps, diff_vec);
CHKERR mField.getInterface<VecManager>()->setOtherLocalGhostVector(
"ADJOINT", "ADJOINT_FIELD", "GEOMETRY", RowColData::ROW,
current_geometry, INSERT_VALUES, SCATTER_REVERSE);
MoFEMFunctionReturn(0);
};
auto fe = get_fd_direvative_fe();
auto fp = vectorDuplicate(diff_vec);
auto fm = vectorDuplicate(diff_vec);
auto calc_impl = [&](auto f, auto eps) { // is this finite difference!
MoFEMFunctionBegin;
CHKERR VecZeroEntries(f);
fe->f = f;
CHKERR perturb_geometry(eps, diff_vec);
CHKERR DMoFEMLoopFiniteElements(simple->getDM(),
simple->getDomainFEName(), fe);
CHKERR VecAssemblyBegin(f);
CHKERR VecAssemblyEnd(f);
CHKERR VecGhostUpdateBegin(f, ADD_VALUES, SCATTER_REVERSE);
CHKERR VecGhostUpdateEnd(f, ADD_VALUES, SCATTER_REVERSE);
auto post_proc_rhs = get_essential_fe();
post_proc_rhs->f = f;
CHKERR DMoFEMPostProcessFiniteElements(simple->getDM(),
post_proc_rhs.get());
MoFEMFunctionReturn(0);
};
CHKERR calc_impl(fp, eps);
CHKERR calc_impl(fm, -eps);
CHKERR VecWAXPY(fd_vec, -1.0, fm, fp);
CHKERR VecScale(fd_vec, 1.0 / (2.0 * eps));
CHKERR mField.getInterface<VecManager>()->setOtherLocalGhostVector(
"ADJOINT", "ADJOINT_FIELD", "GEOMETRY", RowColData::ROW,
initial_current_geometry, INSERT_VALUES, SCATTER_REVERSE);
if constexpr (debug)
CHKERR geom_norm(mField);
MoFEMFunctionReturn(0);
};
// here starts the preperation for the derivative
auto get_direvative_fe = [&](auto diff_vec) {
auto fe_adjoint = boost::make_shared<DomainEle>(mField);
fe_adjoint->getRuleHook = integrationRule;
auto &pip = fe_adjoint->getOpPtrVector();
auto jac_ptr = boost::make_shared<MatrixDouble>();
auto det_ptr = boost::make_shared<VectorDouble>();
auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
auto diff_jac_ptr = boost::make_shared<MatrixDouble>();
auto cof_ptr = boost::make_shared<VectorDouble>();
using OpCoFactor =
AdJoint<DomainEleOp>::Integration<GAUSS>::OpGetCoFactor<SPACE_DIM>;
AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {H1}, "GEOMETRY");
pip.push_back(new OpCalculateVectorFieldGradient<SPACE_DIM, SPACE_DIM>(
"GEOMETRY", jac_ptr));
pip.push_back(new OpCalculateVectorFieldGradient<SPACE_DIM, SPACE_DIM>(
"U", diff_jac_ptr, diff_vec));
pip.push_back(new OpCoFactor(jac_ptr, diff_jac_ptr, cof_ptr));
// Add RHS operators for internal forces
auto common_ptr = HookeOps::commonDataFactory<SPACE_DIM, I, DomainEleOp>(
mField, pip, "U", "MAT_ELASTIC", Sev::noisy);
pip.push_back(new OpAdJointGradTimesSymTensor<SPACE_DIM, I, DomainBaseOp>(
"U", common_ptr, jac_ptr, diff_jac_ptr, cof_ptr));
return fe_adjoint;
};
auto get_objective_fe = [&](auto diff_vec, auto grad_vec,
auto glob_objective_ptr,
auto glob_objective_grad_ptr) {
auto fe_adjoint = boost::make_shared<DomainEle>(mField);
fe_adjoint->getRuleHook = integrationRule;
auto &pip = fe_adjoint->getOpPtrVector();
AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {H1}, "GEOMETRY");
auto jac_ptr = boost::make_shared<MatrixDouble>();
auto det_ptr = boost::make_shared<VectorDouble>();
auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
auto diff_jac_ptr = boost::make_shared<MatrixDouble>();
auto cof_ptr = boost::make_shared<VectorDouble>();
auto d_grad_ptr = boost::make_shared<MatrixDouble>();
auto d_u_ptr = boost::make_shared<MatrixDouble>();
auto u_ptr = boost::make_shared<MatrixDouble>();
using OpCoFactor =
AdJoint<DomainEleOp>::Integration<GAUSS>::OpGetCoFactor<SPACE_DIM>;
pip.push_back(new OpCalculateVectorFieldGradient<SPACE_DIM, SPACE_DIM>(
"GEOMETRY", jac_ptr));
pip.push_back(new OpCalculateVectorFieldGradient<SPACE_DIM, SPACE_DIM>(
"U", diff_jac_ptr,
diff_vec)); // Note: that vector is stored on displacemnt vector, that
// why is used here
pip.push_back(new OpCoFactor(jac_ptr, diff_jac_ptr, cof_ptr));
pip.push_back(new OpCalculateVectorFieldGradient<SPACE_DIM, SPACE_DIM>(
"U", d_grad_ptr, grad_vec));
pip.push_back(
new OpCalculateVectorFieldValues<SPACE_DIM>("U", d_u_ptr, grad_vec));
pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_ptr));
auto common_ptr = HookeOps::commonDataFactory<SPACE_DIM, I, DomainEleOp>(
mField, pip, "U", "MAT_ELASTIC", Sev::noisy);
pip.push_back(new OpAdJointObjective(
pythonPtr, common_ptr, jac_ptr, diff_jac_ptr, cof_ptr, d_grad_ptr,
d_u_ptr, u_ptr, glob_objective_ptr, glob_objective_grad_ptr));
return fe_adjoint;
};
auto dm = simple->getDM();
auto f = createDMVector(dm);
auto d = vectorDuplicate(f);
auto dm_diff_vec = vectorDuplicate(d);
auto zero_diff_vec = vectorDuplicate(dm_diff_vec);
auto zero_state_sensitivity = vectorDuplicate(d);
CHKERR VecZeroEntries(zero_diff_vec);
CHKERR VecZeroEntries(zero_state_sensitivity);
CHKERR VecGhostUpdateBegin(zero_diff_vec, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(zero_diff_vec, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateBegin(zero_state_sensitivity, INSERT_VALUES,
SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(zero_state_sensitivity, INSERT_VALUES,
SCATTER_FORWARD);
auto adjoint_fe = get_direvative_fe(dm_diff_vec);
auto objective_ptr_direct = boost::make_shared<double>(0.0);
auto objective_grad_ptr_direct = boost::make_shared<double>(0.0);
auto objective_fe_direct = get_objective_fe(
dm_diff_vec, d, objective_ptr_direct, objective_grad_ptr_direct);
auto objective_ptr_explicit = boost::make_shared<double>(0.0);
auto objective_grad_ptr_explicit = boost::make_shared<double>(0.0);
auto objective_fe_explicit =
get_objective_fe(dm_diff_vec, zero_state_sensitivity,
objective_ptr_explicit, objective_grad_ptr_explicit);
auto objective_ptr_value = boost::make_shared<double>(0.0);
auto objective_grad_ptr_value = boost::make_shared<double>(0.0);
auto objective_fe_value =
get_objective_fe(zero_diff_vec, zero_state_sensitivity,
objective_ptr_value, objective_grad_ptr_value);
auto set_variance_of_geometry =
[&](auto mode, auto mod_vec) { // think of mod_vec as X(tau) = X + tau*v_h
// take the mod_vec and write it into the adjoint DM then copy that
// field to dm_diff_vec using the same layout as U
MoFEMFunctionBegin;
CHKERR DMoFEMMeshToLocalVector(adjointDM, mod_vec, INSERT_VALUES,
SCATTER_REVERSE);
CHKERR mField.getInterface<VecManager>()->setOtherLocalGhostVector(
simple->getProblemName(), "U", "ADJOINT_FIELD", RowColData::ROW,
dm_diff_vec, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateBegin(dm_diff_vec, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(dm_diff_vec, INSERT_VALUES, SCATTER_FORWARD);
MoFEMFunctionReturn(0);
};
auto calculate_variance_internal_forces = [&](auto mode, auto mod_vec) {
MoFEMFunctionBegin;
CHKERR VecZeroEntries(f);
CHKERR VecGhostUpdateBegin(f, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(f, INSERT_VALUES, SCATTER_FORWARD);
adjoint_fe->f = f;
CHKERR DMoFEMLoopFiniteElements(dm, simple->getDomainFEName(), adjoint_fe);
CHKERR VecAssemblyBegin(f);
CHKERR VecAssemblyEnd(f);
CHKERR VecGhostUpdateBegin(f, ADD_VALUES, SCATTER_REVERSE);
CHKERR VecGhostUpdateEnd(f, ADD_VALUES, SCATTER_REVERSE);
auto post_proc_rhs = get_essential_fe();
post_proc_rhs->f = f;
CHKERR DMoFEMPostProcessFiniteElements(simple->getDM(),
post_proc_rhs.get());
CHKERR VecScale(f, -1.0);
#ifndef NDEBUG
constexpr bool debug = true;
if constexpr (debug) {
double norm0;
CHKERR VecNorm(f, NORM_2, &norm0);
auto fd_check = vectorDuplicate(f);
double eps = 1e-5;
CHKERR calulate_fd_residual(eps, dm_diff_vec, fd_check);
double nrm;
CHKERR VecAXPY(fd_check, -1.0, f);
CHKERR VecNorm(fd_check, NORM_2, &nrm);
MOFEM_LOG("WORLD", Sev::inform)
<< " FD check for internal forces [ " << mode << " ]: " << nrm
<< " / " << norm0 << " ( " << (nrm / norm0) << " )";
}
#endif
MoFEMFunctionReturn(0);
};
auto calculate_variance_of_displacement = [&](auto mode, auto mod_vec) {
MoFEMFunctionBegin;
CHKERR KSPSolve(kspElastic, f, d);
CHKERR VecGhostUpdateBegin(d, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(d, INSERT_VALUES, SCATTER_FORWARD);
MoFEMFunctionReturn(0);
};
auto evaluate_objective_terms =
[&](auto objective_fe, auto objective_ptr, auto objective_grad_ptr,
double &objective_value, double &objective_gradient) {
MoFEMFunctionBegin;
*objective_ptr = 0.0;
*objective_grad_ptr = 0.0;
CHKERR DMoFEMLoopFiniteElements(dm, simple->getDomainFEName(),
objective_fe);
std::array<double, 2> array = {*objective_ptr, *objective_grad_ptr};
MPI_Allreduce(MPI_IN_PLACE, array.data(), 2, MPI_DOUBLE, MPI_SUM,
mField.get_comm());
objective_value = array[0];
objective_gradient = array[1];
MoFEMFunctionReturn(0);
};
auto calculate_objective_value = [&]() {
MoFEMFunctionBegin;
double objective_value = 0;
double objective_gradient = 0;
CHKERR evaluate_objective_terms(objective_fe_value, objective_ptr_value,
objective_grad_ptr_value, objective_value,
objective_gradient);
*objective_function_value = objective_value;
MoFEMFunctionReturn(0);
};
auto calculate_variance_of_objective_function_dJ_du = [&](Vec dJ_du) {
MoFEMFunctionBegin;
auto fe = boost::make_shared<DomainEle>(mField);
fe->getRuleHook = integrationRule;
auto &pip = fe->getOpPtrVector();
AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {H1}, "GEOMETRY");
auto u_ptr = boost::make_shared<MatrixDouble>();
pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_ptr));
auto common_ptr = HookeOps::commonDataFactory<SPACE_DIM, I, DomainEleOp>(
mField, pip, "U", "MAT_ELASTIC", Sev::noisy);
pip.push_back(new OpStateSensitivity("U", pythonPtr, common_ptr, u_ptr));
CHKERR VecZeroEntries(dJ_du);
fe->f = dJ_du;
CHKERR DMoFEMLoopFiniteElements(dm, simple->getDomainFEName(), fe);
CHKERR VecAssemblyBegin(dJ_du);
CHKERR VecAssemblyEnd(dJ_du);
auto post_proc_rhs = get_essential_fe();
post_proc_rhs->f = dJ_du;
CHKERR DMoFEMPostProcessFiniteElements(simple->getDM(),
post_proc_rhs.get());
MoFEMFunctionReturn(0);
};
auto calculate_adjoint_lambda = [&](auto lambda, auto dJ_du) {
MoFEMFunctionBegin;
MOFEM_LOG("WORLD", Sev::inform) << "Solving for adjoint variable lambda";
CHKERR VecZeroEntries(lambda);
CHKERR KSPSolveTranspose(kspElastic, dJ_du, lambda);
CHKERR VecGhostUpdateBegin(lambda, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(lambda, INSERT_VALUES, SCATTER_FORWARD);
MoFEMFunctionReturn(0);
};
CHKERR VecZeroEntries(objective_function_gradient);
CHKERR VecZeroEntries(adjoint_vector);
CHKERR calculate_objective_value();
MOFEM_LOG("WORLD", Sev::verbose)
<< "Objective function: " << *objective_function_value;
auto direct = [&]() {
MoFEMFunctionBegin;
int mode = 0;
for (auto mod_vec : modeVecs) {
CHKERR set_variance_of_geometry(mode, mod_vec);
CHKERR calculate_variance_internal_forces(mode, mod_vec);
CHKERR calculate_variance_of_displacement(mode, mod_vec);
double objective_value = 0;
double objective_gradient = 0;
CHKERR evaluate_objective_terms(objective_fe_direct, objective_ptr_direct,
objective_grad_ptr_direct,
objective_value, objective_gradient);
CHKERR VecSetValue(objective_function_gradient, mode, objective_gradient,
INSERT_VALUES);
CHKERR VecAXPY(adjoint_vector, objective_gradient, dm_diff_vec);
++mode;
}
MoFEMFunctionReturn(0);
};
auto lambda = vectorDuplicate(f);
auto dJ_du = vectorDuplicate(f);
auto adjoint = [&]() {
MoFEMFunctionBegin;
CHKERR calculate_variance_of_objective_function_dJ_du(dJ_du);
CHKERR calculate_adjoint_lambda(lambda, dJ_du);
int mode = 0;
for (auto mod_vec : modeVecs) {
CHKERR set_variance_of_geometry(mode, mod_vec);
CHKERR calculate_variance_internal_forces(mode, mod_vec);
double objective_value = 0;
double objective_gradient_explicit = 0;
CHKERR evaluate_objective_terms(
objective_fe_explicit, objective_ptr_explicit,
objective_grad_ptr_explicit, objective_value,
objective_gradient_explicit);
double lambda_dot_residual_variation = 0;
CHKERR VecDot(lambda, f, &lambda_dot_residual_variation);
const double dJ_dp =
objective_gradient_explicit + lambda_dot_residual_variation;
CHKERR VecSetValue(objective_function_gradient, mode, dJ_dp,
INSERT_VALUES);
CHKERR VecAXPY(adjoint_vector, dJ_dp, dm_diff_vec);
++mode;
}
CHKERR VecAssemblyBegin(objective_function_gradient);
CHKERR VecAssemblyEnd(objective_function_gradient);
MoFEMFunctionReturn(0);
};
switch (derivative_type) {
case DIRECT:
MOFEM_LOG("WORLD", Sev::inform) << "Running Direct Sensitivity...";
CHKERR direct();
break;
case ADJOINT:
MOFEM_LOG("WORLD", Sev::inform) << "Running Adjoint Sensitivity...";
CHKERR adjoint();
break;
default:
SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
"Wrong sensitivity type selected");
}
CHKERR VecAssemblyBegin(objective_function_gradient);
CHKERR VecAssemblyEnd(objective_function_gradient);
CHKERR VecAssemblyBegin(adjoint_vector);
CHKERR VecAssemblyEnd(adjoint_vector);
CHKERR VecGhostUpdateBegin(adjoint_vector, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(adjoint_vector, INSERT_VALUES, SCATTER_FORWARD);
MoFEMFunctionReturn(0);
}
//! [calculateGradient]
static char help[] = "...\n\n";
/**
* @brief Main function for topology optimization tutorial using adjoint method
*
* This tutorial demonstrates structural topology optimization using:
* - MoFEM finite element library for elasticity analysis
* - Adjoint method for efficient gradient computation
* - TAO optimization library for design optimization
* - Python interface for flexible objective function definition
*
* Workflow:
* 1. Initialize MoFEM, PETSc and Python environments
* 2. Read mesh and setup finite element problem
* 3. Define objective function via Python interface
* 4. Compute topology optimization modes as design variables
* 5. Run gradient-based optimization using adjoint sensitivities
* 6. Post-process optimized design
*
* The adjoint method enables efficient gradient computation with cost
* independent of the number of design variables, making it suitable
* for large-scale topology optimization problems.
*
* Required input files:
* - Mesh file (.h5m format from CUBIT)
* - Parameter file (param_file.petsc)
* - Objective function (objective_function.py)
*
* @param argc Command line argument count
* @param argv Command line argument values
* @return int Exit code (0 for success)
*/
int main(int argc, char *argv[]) {
// Initialize Python environment for objective function interface
Py_Initialize();
np::initialize();
// Initialize MoFEM/PETSc and MOAB data structures
const char param_file[] = "param_file.petsc";
MoFEM::Core::Initialize(&argc, &argv, param_file, help);
auto core_log = logging::core::get();
core_log->add_sink(
LogManager::createSink(LogManager::getStrmWorld(), "TIMER"));
core_log->add_sink(
LogManager::createSink(LogManager::getStrmSync(), "FieldEvaluator"));
LogManager::setLog("FieldEvaluator");
MOFEM_LOG_TAG("FieldEvaluator", "field_eval");
try {
//! [Register MoFEM discrete manager in PETSc]
DMType dm_name = "DMMOFEM";
CHKERR DMRegister_MoFEM(dm_name);
DMType dm_name_mg = "DMMOFEM_MG";
CHKERR DMRegister_MGViaApproxOrders(dm_name_mg);
//! [Register MoFEM discrete manager in PETSc
//! [Create MoAB]
moab::Core mb_instance; ///< mesh database
moab::Interface &moab = mb_instance; ///< mesh database interface
//! [Create MoAB]
//! [Create MoFEM]
MoFEM::Core core(moab); ///< finite element database
MoFEM::Interface &m_field = core; ///< finite element database interface
//! [Create MoFEM]
//! [Example]
Example ex(m_field);
CHKERR ex.runProblem();
//! [Example]
}
CATCH_ERRORS;
CHKERR MoFEM::Core::Finalize();
if (Py_FinalizeEx() < 0) {
exit(120);
}
}
Range get_range_from_block(MoFEM::Interface &m_field,
const std::string block_name, int dim) {
Range r;
auto mesh_mng = m_field.getInterface<MeshsetsManager>();
auto bcs = mesh_mng->getCubitMeshsetPtr(
std::regex((boost::format("%s(.*)") % block_name).str())
);
for (auto bc : bcs) {
Range faces;
CHK_MOAB_THROW(bc->getMeshsetIdEntitiesByDimension(m_field.get_moab(), dim,
faces, true),
"get meshset ents");
r.merge(faces);
}
for (auto dd = dim - 1; dd >= 0; --dd) {
if (dd >= 0) {
Range ents;
CHK_MOAB_THROW(m_field.get_moab().get_adjacencies(r, dd, false, ents,
moab::Interface::UNION),
"get adjs");
r.merge(ents);
} else {
Range verts;
CHK_MOAB_THROW(m_field.get_moab().get_connectivity(r, verts),
"get verts");
r.merge(verts);
}
CHK_MOAB_THROW(
m_field.getInterface<CommInterface>()->synchroniseEntities(r), "comm");
}
return r;
};
MoFEMErrorCode save_range(moab::Interface &moab, const std::string name,
const Range r) {
MoFEMFunctionBegin;
auto out_meshset = get_temp_meshset_ptr(moab);
CHKERR moab.add_entities(*out_meshset, r);
CHKERR moab.write_file(name.c_str(), "VTK", "", out_meshset->get_ptr(), 1);
MoFEMFunctionReturn(0);
};
CalculateTraction.hpp
Calculate traction for linear problem.
ElasticSpring.hpp
Implementation of elastic spring bc.
FluidLevel.hpp
Natural boundary condition applying pressure from fluid.
HookeOps.hpp
Implementation of Hookes operator Hookes for linear elastic problems in MoFEM.
type
std::string type
Definition JsonConfigManager.cpp:28
MOFEM_LOG_SYNCHRONISE
#define MOFEM_LOG_SYNCHRONISE(comm)
Synchronise "SYNC" channel.
Definition LogManager.hpp:353
NaturalBoundaryBC.hpp
Implementation of natural boundary conditions.
NaturalDomainBC.hpp
Boundary conditions in domain, i.e. body forces.
ObjectiveFunctionData.hpp
Interface for Python-based objective function evaluation in topology optimization.
size
PetscInt size
Definition TaoNullSpace.cpp:125
FTENSOR_INDEX
#define FTENSOR_INDEX(DIM, I)
Definition Templates.hpp:2668
name
std::string name
Definition TestingInterface.cpp:88
help
static char help[]
Definition activate_deactivate_dofs.cpp:13
derivative_type
SensitivityMethod derivative_type
Definition adjoint.cpp:109
SensitivityMethod
SensitivityMethod
Definition adjoint.cpp:107
DIRECT
@ DIRECT
Definition adjoint.cpp:107
ADJOINT
@ ADJOINT
Definition adjoint.cpp:107
DomainBaseOp
FormsIntegrators< DomainEleOp >::Assembly< A >::OpBase DomainBaseOp
[Postprocess results]
Definition adjoint.cpp:1383
get_range_from_block
Range get_range_from_block(MoFEM::Interface &m_field, const std::string block_name, int dim)
Definition adjoint.cpp:2288
main
int main()
Definition adol-c_atom.cpp:46
a
constexpr double a
Definition approx_sphere.cpp:30
SPACE_DIM
constexpr int SPACE_DIM
Definition child_and_parent.cpp:17
DomainEle
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
Definition child_and_parent.cpp:35
BoundaryEle
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
Definition child_and_parent.cpp:40
QUIET
@ QUIET
Definition definitions.h:221
ROW
@ ROW
Definition definitions.h:136
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition definitions.h:385
MF_EXIST
@ MF_EXIST
Definition definitions.h:113
FieldApproximationBase
FieldApproximationBase
approximation base
Definition definitions.h:58
LASTBASE
@ LASTBASE
Definition definitions.h:69
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition definitions.h:60
DEMKOWICZ_JACOBI_BASE
@ DEMKOWICZ_JACOBI_BASE
Definition definitions.h:66
CHK_THROW_MESSAGE
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
Definition definitions.h:612
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition definitions.h:463
H1
@ H1
continuous field
Definition definitions.h:85
NOSPACE
@ NOSPACE
Definition definitions.h:83
MYPCOMM_INDEX
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition definitions.h:228
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition definitions.h:359
CHK_MOAB_THROW
#define CHK_MOAB_THROW(err, msg)
Check error code of MoAB function and throw MoFEM exception.
Definition definitions.h:592
MOFEM_NOT_FOUND
@ MOFEM_NOT_FOUND
Definition definitions.h:33
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition definitions.h:31
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition definitions.h:432
CHKERR
#define CHKERR
Inline error check.
Definition definitions.h:551
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition definitions.h:456
BASE_DIM
constexpr int BASE_DIM
Definition dg_projection.cpp:35
order
constexpr int order
Definition dg_projection.cpp:38
debug
static const bool debug
Definition dm_create_subdm.cpp:12
PostProcEleDomain
PostProcEleByDim< SPACE_DIM >::PostProcEleDomain PostProcEleDomain
Definition dynamic_first_order_con_law.cpp:51
PostProcEleBdy
PostProcEleByDim< SPACE_DIM >::PostProcEleBdy PostProcEleBdy
Definition dynamic_first_order_con_law.cpp:53
bulk_modulus_K
double bulk_modulus_K
Definition dynamic_first_order_con_law.cpp:95
shear_modulus_G
double shear_modulus_G
Definition dynamic_first_order_con_law.cpp:96
F
@ F
Definition free_surface.cpp:630
is_plane_strain
PetscBool is_plane_strain
Definition gradient.cpp:37
DomainEle
PipelineManager::ElementsAndOpsByDim< SPACE_DIM >::DomainEle DomainEle
Domain finite elements.
Definition gradient.cpp:45
SPACE_DIM
constexpr int SPACE_DIM
[Define dimension]
Definition gradient.cpp:19
derivative_type
SensitivityMethod derivative_type
Definition gradient.cpp:96
I
constexpr IntegrationType I
Use Gauss quadrature for integration.
Definition gradient.cpp:25
DIRECT
@ DIRECT
Definition gradient.cpp:94
ADJOINT
@ ADJOINT
Definition gradient.cpp:94
DomainBaseOp
FormsIntegrators< DomainEleOp >::Assembly< A >::OpBase DomainBaseOp
[Postprocess results]
Definition gradient.cpp:699
A
constexpr AssemblyType A
[Define dimension]
Definition gradient.cpp:23
BoundaryEle
PipelineManager::ElementsAndOpsByDim< SPACE_DIM >::BoundaryEle BoundaryEle
Boundary finite elements.
Definition gradient.cpp:47
get_range_from_block
Range get_range_from_block(MoFEM::Interface &m_field, const std::string block_name, int dim)
Definition gradient.cpp:1347
MoFEM::DMMoFEMSetIsPartitioned
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
Definition DMMoFEM.cpp:1142
MoFEM::DMMoFEMCreateSubDM
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:244
MoFEM::DMMoFEMAddElement
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
Definition DMMoFEM.cpp:517
MoFEM::DMMoFEMSetSquareProblem
PetscErrorCode DMMoFEMSetSquareProblem(DM dm, PetscBool square_problem)
set squared problem
Definition DMMoFEM.cpp:479
MoFEM::DMMoFEMCreateMoFEM
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:143
MoFEM::DMoFEMPostProcessFiniteElements
PetscErrorCode DMoFEMPostProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
Definition DMMoFEM.cpp:575
MoFEM::DMMoFEMAddSubFieldRow
PetscErrorCode DMMoFEMAddSubFieldRow(DM dm, const char field_name[])
Definition DMMoFEM.cpp:267
MoFEM::DMoFEMMeshToLocalVector
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode, RowColData rc=RowColData::COL)
set local (or ghosted) vector values on mesh for partition only
Definition DMMoFEM.cpp:543
MoFEM::DMRegister_MoFEM
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition DMMoFEM.cpp:43
MoFEM::DMRegister_MGViaApproxOrders
MoFEMErrorCode DMRegister_MGViaApproxOrders(const char sname[])
Register DM for Multi-Grid via approximation orders.
Definition PCMGSetUpViaApproxOrders.cpp:302
MoFEM::DMoFEMLoopFiniteElements
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:605
MoFEM::createDMVector
auto createDMVector(DM dm, RowColData rc=RowColData::COL)
Get smart vector from DM.
Definition DMMoFEM.hpp:1250
MoFEM::DMMoFEMAddSubFieldCol
PetscErrorCode DMMoFEMAddSubFieldCol(DM dm, const char field_name[])
Definition DMMoFEM.cpp:309
MoFEM::createDMMatrix
auto createDMMatrix(DM dm)
Get smart matrix from DM.
Definition DMMoFEM.hpp:1207
MoFEM::DMoFEMPreProcessFiniteElements
PetscErrorCode DMoFEMPreProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
Definition DMMoFEM.cpp:565
MoFEM::CoreInterface::add_ents_to_finite_element_by_dim
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
MoFEM::CoreInterface::add_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
MoFEM::CoreInterface::build_finite_elements
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
MoFEM::CoreInterface::modify_finite_element_add_field_col
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
MoFEM::CoreInterface::modify_finite_element_add_field_row
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
MoFEM::CoreInterface::modify_finite_element_add_field_data
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string name_field)=0
set finite element field data
MoFEM::CoreInterface::build_fields
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
MoFEM::CoreInterface::add_ents_to_field_by_dim
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.
MoFEM::CoreInterface::set_field_order
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.
MoFEM::IntegrationType
IntegrationType
Form integrator integration types.
Definition FormsIntegrators.hpp:237
MoFEM::AssemblyType
AssemblyType
[Storage and set boundary conditions]
Definition FormsIntegrators.hpp:200
MoFEM::GAUSS
@ GAUSS
Gaussian quadrature integration.
Definition FormsIntegrators.hpp:238
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition LogManager.hpp:316
MOFEM_LOG_TAG
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition LogManager.hpp:347
MOFEM_LOG_CHANNEL
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
Definition LogManager.hpp:292
MoFEM::CoreInterface::loop_dofs
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.
MoFEM::MeshsetsManager::getCubitMeshsetPtr
MoFEMErrorCode getCubitMeshsetPtr(const int ms_id, const CubitBCType cubit_bc_type, const CubitMeshSets **cubit_meshset_ptr) const
get cubit meshset
Definition MeshsetsManager.cpp:745
i
FTensor::Index< 'i', SPACE_DIM > i
Definition hcurl_divergence_operator_2d.cpp:27
lambda
static double lambda
Definition incompressible_elasticity.cpp:181
testing
PetscBool testing
Definition initial_diffusion.cpp:62
c
const double c
speed of light (cm/ns)
Definition initial_diffusion.cpp:40
v
const double v
phase velocity of light in medium (cm/ns)
Definition initial_diffusion.cpp:41
J
FTensor::Index< 'J', DIM1 > J
Definition level_set.cpp:30
l
FTensor::Index< 'l', 3 > l
Definition matrix_function.cpp:23
j
FTensor::Index< 'j', 3 > j
Definition matrix_function.cpp:21
k
FTensor::Index< 'k', 3 > k
Definition matrix_function.cpp:22
EigenMatrix::Vec
const FTensor::Tensor2< T, Dim, Dim > Vec
Definition MatrixFunction.hpp:64
ElasticPostProc::postProcessElasticResults
MoFEMErrorCode postProcessElasticResults(MoFEM::Interface &mField, SmartPetscObj< DM > dm, const std::string &domain_fe_name, const std::string &out_file_name, std::vector< std::pair< std::string, SmartPetscObj< Vec > > > extra_vectors={}, const std::vector< std::string > &tags_to_transfer={}, const Sev hooke_ops_sev=Sev::verbose)
Definition ElasticPostProc.hpp:10
FTensor::dd
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition ddTensor0.hpp:33
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition Exceptions.hpp:56
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition Types.hpp:77
MoFEM::Types::VectorDouble
UBlasVector< double > VectorDouble
Definition Types.hpp:68
MoFEM
implementation of Data Operators for Forces and Sources
Definition Common.hpp:10
MoFEM::checkValue< TestType::DERIVATIVE >
MoFEMErrorCode TestingInterface::checkValue< TestType::DERIVATIVE >(const std::string &check_name, PetscReal actual, PetscReal expected, Tolerance tolerance, Sev error_severity, Sev success_severity) const
Definition TestingInterface.cpp:747
MoFEM::createKSP
auto createKSP(MPI_Comm comm)
Definition PetscSmartObj.hpp:269
MoFEM::DMMoFEMSetDestroyProblem
PetscErrorCode DMMoFEMSetDestroyProblem(DM dm, PetscBool destroy_problem)
Definition DMMoFEM.cpp:463
MoFEM::PetscOptionsGetInt
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
Definition DeprecatedPetsc.hpp:142
MoFEM::PetscOptionsGetReal
PetscErrorCode PetscOptionsGetReal(PetscOptions *, const char pre[], const char name[], PetscReal *dval, PetscBool *set)
Definition DeprecatedPetsc.hpp:152
MoFEM::PetscOptionsGetBool
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
Definition DeprecatedPetsc.hpp:182
MoFEM::vectorDuplicate
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
Definition PetscSmartObj.hpp:223
MoFEM::PetscOptionsGetRealArray
PetscErrorCode PetscOptionsGetRealArray(PetscOptions *, const char pre[], const char name[], PetscReal dval[], PetscInt *nmax, PetscBool *set)
Definition DeprecatedPetsc.hpp:192
MoFEM::createVectorMPI
auto createVectorMPI(MPI_Comm comm, PetscInt n, PetscInt N)
Create MPI Vector.
Definition PetscSmartObj.hpp:204
MoFEM::getDMKspCtx
auto getDMKspCtx(DM dm)
Get KSP context data structure used by DM.
Definition DMMoFEM.hpp:1264
MoFEM::invertTensor
static MoFEMErrorCode invertTensor(FTensor::Tensor2< T1, DIM, DIM > &t, T2 &det, FTensor::Tensor2< T3, DIM, DIM > &inv_t)
Definition Templates.hpp:2411
MoFEM::determinantTensor
static auto determinantTensor(FTensor::Tensor2< T, DIM, DIM > &t)
Calculate the determinant of a tensor of rank DIM.
Definition Templates.hpp:2254
MoFEM::PetscOptionsGetEList
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
Definition DeprecatedPetsc.hpp:203
MoFEM::createGhostVector
auto createGhostVector(MPI_Comm comm, PetscInt n, PetscInt N, PetscInt nghost, const PetscInt ghosts[])
Create smart ghost vector.
Definition PetscSmartObj.hpp:181
MoFEM::PetscOptionsGetString
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
Definition DeprecatedPetsc.hpp:172
MoFEM::getFTensor0FromMat
auto getFTensor0FromMat(M &data)
Get tensor rank 0 (scalar) form data vector.
Definition Templates.hpp:197
MoFEM::get_temp_meshset_ptr
auto get_temp_meshset_ptr(moab::Interface &moab)
Create smart pointer to temporary meshset.
Definition Templates.hpp:2531
MoFEM::TaoSetObjectiveAndGradient
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:386
MoFEM::getFTensor0FromVec
static auto getFTensor0FromVec(V &data)
Get tensor rank 0 (scalar) form data vector.
Definition Templates.hpp:163
MoFEM::createDM
auto createDM(MPI_Comm comm, const std::string dm_type_name)
Creates smart DM object.
Definition PetscSmartObj.hpp:143
MoFEM::VecSetValues
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
Definition EntitiesFieldData.hpp:1874
MoFEM::createTao
auto createTao(MPI_Comm comm)
Definition PetscSmartObj.hpp:251
ShapeOptimization::create_python_objective_function
boost::shared_ptr< ObjectiveFunctionData > create_python_objective_function(std::string)
Definition ObjectiveFunctionData.cpp:1322
profile_test.file
file
Definition profile_test.py:2163
sdf_hertz_2d_axisymm_plane.d
float d
Definition sdf_hertz_2d_axisymm_plane.py:4
sdf.r
int r
Definition sdf.py:205
I
constexpr IntegrationType I
Definition operators_tests.cpp:31
A
constexpr AssemblyType A
Definition operators_tests.cpp:30
OpPPMap
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
Definition photon_diffusion.cpp:40
young_modulus
double young_modulus
Young modulus.
Definition plastic.cpp:126
BoundaryLhsBCs
NaturalBC< BoundaryEleOp >::Assembly< AT >::BiLinearForm< IT > BoundaryLhsBCs
Definition plastic.cpp:176
BoundaryRhsBCs
NaturalBC< BoundaryEleOp >::Assembly< AT >::LinearForm< IT > BoundaryRhsBCs
Definition plastic.cpp:173
EXECUTABLE_DIMENSION
#define EXECUTABLE_DIMENSION
Definition plastic.cpp:13
do_eval_field
PetscBool do_eval_field
Evaluate field.
Definition plastic.cpp:120
OpBoundaryLhsBCs
BoundaryLhsBCs::OpFlux< PlasticOps::BoundaryBCs, 1, SPACE_DIM > OpBoundaryLhsBCs
Definition plastic.cpp:178
poisson_ratio
double poisson_ratio
Poisson ratio.
Definition plastic.cpp:127
DomainRhsBCs
NaturalBC< DomainEleOp >::Assembly< AT >::LinearForm< IT > DomainRhsBCs
Definition plastic.cpp:170
OpDomainRhsBCs
DomainRhsBCs::OpFlux< PlasticOps::DomainBCs, 1, SPACE_DIM > OpDomainRhsBCs
Definition plastic.cpp:172
OpBoundaryRhsBCs
BoundaryRhsBCs::OpFlux< PlasticOps::BoundaryBCs, 1, SPACE_DIM > OpBoundaryRhsBCs
Definition plastic.cpp:175
SideEle
ElementsAndOps< SPACE_DIM >::SideEle SideEle
Definition plastic.cpp:62
field_name
constexpr auto field_name
Definition poisson_2d_homogeneous.cpp:13
OpBase
OpBaseImpl< PETSC, EdgeEleOp > OpBase
Definition radiation.cpp:29
OpMass
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
[Only used with Hooke equation (linear material model)]
Definition seepage.cpp:56
is_plane_strain
PetscBool is_plane_strain
Definition seepage.cpp:177
g
constexpr double g
Definition shallow_wave.cpp:63
m
FTensor::Index< 'm', 3 > m
Definition shallow_wave.cpp:80
BoundaryBCs
Boundary conditions marker.
Definition arc_length_interface.cpp:23
DomainBCs
[Define entities]
Definition elastic.cpp:38
Example
[Example]
Definition plastic.cpp:217
Example::boundaryCondition
MoFEMErrorCode boundaryCondition()
[Set up problem]
Definition dynamic_first_order_con_law.cpp:532
Example::assembleSystem
MoFEMErrorCode assembleSystem()
[Push operators to pipeline]
Definition dynamic_first_order_con_law.cpp:635
Example::readMesh
MoFEMErrorCode readMesh()
[Run problem]
Definition dynamic_first_order_con_law.cpp:459
Example::integrationRule
static int integrationRule(int, int, int p_data)
Definition radiation.cpp:52
Example::pythonPtr
boost::shared_ptr< ObjectiveFunctionData > pythonPtr
Interface to Python objective function.
Definition adjoint.cpp:170
Example::modeVecs
std::vector< SmartPetscObj< Vec > > modeVecs
Topology mode vectors (design variables)
Definition adjoint.cpp:174
Example::adjointDM
SmartPetscObj< DM > adjointDM
Data manager for adjoint problem.
Definition adjoint.cpp:168
Example::base
FieldApproximationBase base
Choice of finite element basis functions.
Definition plot_base.cpp:68
Example::modeCentroids
std::vector< std::array< double, 3 > > modeCentroids
Centroids of optimization blocks.
Definition adjoint.cpp:176
Example::topologyModes
MoFEMErrorCode topologyModes()
Compute topology optimization modes.
Definition adjoint.cpp:775
Example::simple
Simple * simple
Definition schrod_eig.cpp:63
Example::kspElastic
SmartPetscObj< KSP > kspElastic
Linear solver for elastic problem.
Definition adjoint.cpp:167
Example::M
SmartPetscObj< Mat > M
Definition schrod_eig.cpp:73
Example::initialGeometry
SmartPetscObj< Vec > initialGeometry
Initial geometry field.
Definition adjoint.cpp:179
Example::fieldOrder
int fieldOrder
Polynomial order for approximation.
Definition adjoint.cpp:164
Example::runProblem
MoFEMErrorCode runProblem()
[Run problem]
Definition plastic.cpp:255
Example::calculateGradient
MoFEMErrorCode calculateGradient(PetscReal *objective_function_value, Vec objective_function_gradient, Vec adjoint_vector)
Calculate objective function gradient using adjoint method.
Definition adjoint.cpp:1768
Example::setupAdJoint
MoFEMErrorCode setupAdJoint()
Setup adjoint fields and finite elements.
Definition adjoint.cpp:652
Example::mField
MoFEM::Interface & mField
Reference to MoFEM interface.
Definition plastic.cpp:227
Example::modeBBoxes
std::vector< std::array< double, 6 > > modeBBoxes
Bounding boxes of optimization blocks.
Definition adjoint.cpp:178
Example::setupProblem
MoFEMErrorCode setupProblem()
[Run problem]
Definition plastic.cpp:274
Example::postprocessElastic
MoFEMErrorCode postprocessElastic(int iter, SmartPetscObj< Vec > adjoint_vector=nullptr)
Post-process and output results.
Definition adjoint.cpp:1368
Example::eps
SmartPetscObj< EPS > eps
Definition schrod_eig.cpp:75
Example::solveElastic
MoFEMErrorCode solveElastic()
Solve forward elastic problem.
Definition adjoint.cpp:1240
MoFEM::AddHOOps
Add operators pushing bases from local to physical configuration.
Definition HODataOperators.hpp:642
MoFEM::BcManager
Boundary condition manager for finite element problem setup.
Definition BcManager.hpp:384
MoFEM::CommInterface
Managing BitRefLevels.
Definition CommInterface.hpp:66
MoFEM::CommInterface::synchroniseEntities
MoFEMErrorCode synchroniseEntities(Range &ent, std::map< int, Range > *received_ents, int verb=DEFAULT_VERBOSITY)
synchronize entity range on processors (collective)
Definition CommInterface.cpp:19
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
MoFEM::CoreInterface::build_adjacencies
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
MoFEM::CoreInterface::add_field
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.
MoFEM::CoreInterface::get_comm
virtual MPI_Comm & get_comm() const =0
MoFEM::CoreInterface::get_comm_rank
virtual int get_comm_rank() const =0
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition Core.hpp:83
MoFEM::CoreTmp< 0 >::Initialize
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition Core.cpp:68
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition Core.cpp:123
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition DeprecatedCoreInterface.hpp:16
MoFEM::DisplacementCubitBcData
Definition of the displacement bc data structure.
Definition BCData.hpp:72
MoFEM::EntitiesFieldData::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition EntitiesFieldData.hpp:150
MoFEM::EntitiesFieldData::EntData::getFTensor0N
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
Definition EntitiesFieldData.hpp:1806
MoFEM::EntitiesFieldData::EntData::getFTensor1DiffN
auto getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
Definition EntitiesFieldData.hpp:802
MoFEM::EntitiesFieldData::EntData::getN
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
Definition EntitiesFieldData.hpp:1600
MoFEM::EntitiesFieldData::EntData::getIndices
const VectorInt & getIndices() const
Get global indices of degrees of freedom on entity.
Definition EntitiesFieldData.hpp:1523
MoFEM::EssentialPostProcLhs
Class (Function) to enforce essential constrains on the left hand side diagonal.
Definition Essential.hpp:33
MoFEM::EssentialPostProcRhs
Class (Function) to enforce essential constrains on the right hand side diagonal.
Definition Essential.hpp:41
MoFEM::EssentialPreProc
Class (Function) to enforce essential constrains.
Definition Essential.hpp:25
MoFEM::FieldBlas
Basic algebra on fields.
Definition FieldBlas.hpp:21
MoFEM::FieldBlas::fieldLambdaOnValues
MoFEMErrorCode fieldLambdaOnValues(OneFieldFunctionOnValues lambda, const std::string field_name, Range *ents_ptr=nullptr)
field lambda
Definition FieldBlas.cpp:21
MoFEM::FieldEvaluatorInterface
Field evaluator interface.
Definition FieldEvaluator.hpp:21
MoFEM::MeshsetsManager
Interface for managing meshsets containing materials and boundary conditions.
Definition MeshsetsManager.hpp:208
MoFEM::NaturalBC::Assembly
Assembly methods.
Definition Natural.hpp:65
MoFEM::OpBaseImpl::nbRows
int nbRows
number of dofs on rows
Definition FormsIntegrators.hpp:369
MoFEM::OpBaseImpl::nbRowBaseFunctions
int nbRowBaseFunctions
number or row base functions
Definition FormsIntegrators.hpp:372
MoFEM::OpCalculateVectorFieldGradient
Get field gradients at integration pts for scalar field rank 0, i.e. vector field.
Definition UserDataOperators.hpp:2008
MoFEM::OpCalculateVectorFieldValues
Specialization for MatrixDouble vector field values calculation.
Definition UserDataOperators.hpp:685
MoFEM::OpPostProcMapInMoab
Post post-proc data at points from hash maps.
Definition PostProcBrokenMeshInMoabBase.hpp:745
MoFEM::OperatorsTester
Calculate directional derivative of the right hand side and compare it with tangent matrix derivative...
Definition OperatorsTester.hpp:21
MoFEM::PipelineManager::ElementsAndOpsByDim
Template struct for dimension-specific finite element types.
Definition PipelineManager.hpp:55
MoFEM::PipelineManager
PipelineManager interface.
Definition PipelineManager.hpp:24
MoFEM::PipelineManager::setDomainRhsIntegrationRule
MoFEMErrorCode setDomainRhsIntegrationRule(RuleHookFun rule)
Set integration rule for domain right-hand side finite element.
Definition PipelineManager.hpp:1019
MoFEM::ProblemsManager
Problem manager is used to build and partition problems.
Definition ProblemsManager.hpp:133
MoFEM::Projection10NodeCoordsOnField
Projection of edge entities with one mid-node on hierarchical basis.
Definition Projection10NodeCoordsOnField.hpp:24
MoFEM::Simple
Simple interface for fast problem set-up.
Definition Simple.hpp:27
MoFEM::Simple::addDomainField
MoFEMErrorCode addDomainField(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_ZERO, int verb=-1)
Add field on domain.
Definition Simple.cpp:261
MoFEM::Simple::getMeshset
EntityHandle & getMeshset()
Get the MeshSet object.
Definition Simple.hpp:394
MoFEM::Simple::loadFile
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name, LoadFileFunc loadFunc=defaultLoadFileFunc)
Load mesh file.
Definition Simple.cpp:191
MoFEM::Simple::addBoundaryField
MoFEMErrorCode addBoundaryField(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_ZERO, int verb=-1)
Add field on boundary.
Definition Simple.cpp:355
MoFEM::Simple::getOptions
MoFEMErrorCode getOptions()
get options
Definition Simple.cpp:180
MoFEM::Simple::getDM
MoFEMErrorCode getDM(DM *dm)
Get DM.
Definition Simple.cpp:799
MoFEM::Simple::setFieldOrder
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
Definition Simple.cpp:575
MoFEM::Simple::getBitRefLevelMask
BitRefLevel & getBitRefLevelMask()
Get the BitRefLevelMask.
Definition Simple.hpp:422
MoFEM::Simple::addDataField
MoFEMErrorCode addDataField(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_ZERO, int verb=-1)
Add data field.
Definition Simple.cpp:393
MoFEM::Simple::setUp
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
Definition Simple.cpp:735
MoFEM::Simple::getProblemName
const std::string getProblemName() const
Get the Problem Name.
Definition Simple.hpp:450
MoFEM::Simple::getDomainFEName
const std::string getDomainFEName() const
Get the Domain FE Name.
Definition Simple.hpp:429
MoFEM::Simple::getBitRefLevel
BitRefLevel & getBitRefLevel()
Get the BitRefLevel.
Definition Simple.hpp:415
MoFEM::SmartPetscObj
intrusive_ptr for managing petsc objects
Definition PetscSmartObj.hpp:85
MoFEM::TestingInterface::addTestCommandLine
virtual MoFEMErrorCode addTestCommandLine(const std::string &prefix, int number, const std::string &name, TestHook hook)=0
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition UnknownInterface.hpp:103
MoFEM::VecManager
Vector manager is used to create vectors \mofem_vectors.
Definition VecManager.hpp:23
OpAdJointGradTimesSymTensor
Forward declaration of operator for gradient times symmetric tensor operations.
Definition adjoint.cpp:105
OpAdJointObjective::dGradPtr
boost::shared_ptr< MatrixDouble > dGradPtr
Definition adjoint.cpp:1760
OpAdJointObjective::globObjectiveGradPtr
boost::shared_ptr< double > globObjectiveGradPtr
Definition adjoint.cpp:1765
OpAdJointObjective::commPtr
boost::shared_ptr< HookeOps::CommonData > commPtr
Definition adjoint.cpp:1756
OpAdJointObjective::globObjectivePtr
boost::shared_ptr< double > globObjectivePtr
Definition adjoint.cpp:1764
OpAdJointObjective::pythonPtr
boost::shared_ptr< ObjectiveFunctionData > pythonPtr
Definition adjoint.cpp:1755
OpAdJointObjective::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Compute objective function contributions at element level.
Definition adjoint.cpp:1618
OpAdJointObjective::cofVals
boost::shared_ptr< VectorDouble > cofVals
Definition adjoint.cpp:1759
OpAdJointObjective::uPtr
boost::shared_ptr< MatrixDouble > uPtr
Definition adjoint.cpp:1762
OpAdJointObjective::jacPtr
boost::shared_ptr< MatrixDouble > jacPtr
Definition adjoint.cpp:1757
OpAdJointObjective::diffJacPtr
boost::shared_ptr< MatrixDouble > diffJacPtr
Definition adjoint.cpp:1758
OpAdJointObjective::dUPtr
boost::shared_ptr< MatrixDouble > dUPtr
Definition adjoint.cpp:1761
OpStateSensitivity::commPtr
boost::shared_ptr< HookeOps::CommonData > commPtr
Definition adjoint.cpp:1590
OpStateSensitivity::iNtegrate
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data)
Class dedicated to integrate operator.
Definition adjoint.cpp:1512
OpStateSensitivity::uPtr
boost::shared_ptr< MatrixDouble > uPtr
Definition adjoint.cpp:1591
OpStateSensitivity::pythonPtr
boost::shared_ptr< ObjectiveFunctionData > pythonPtr
Definition adjoint.cpp:1589
save_range
auto save_range
Definition thermo_elastic.cpp:111
Generated by Doxygen 1.9.8 and hosted at University of Glasgow