diff --git a/python/py_akantu.cc b/python/py_akantu.cc
index 071814a72..245d02bd6 100644
--- a/python/py_akantu.cc
+++ b/python/py_akantu.cc
@@ -1,169 +1,169 @@
/**
* Copyright (©) 2018-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* This file is part of Akantu
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see .
*/
/* -------------------------------------------------------------------------- */
#include "aka_config.hh"
// for NLSNotConvergedException
#include "non_linear_solver.hh"
/* -------------------------------------------------------------------------- */
#include "py_aka_common.hh"
#include "py_aka_error.hh"
#include "py_boundary_conditions.hh"
#include "py_constitutive_law.hh"
#include "py_constitutive_law_selector.hh"
#include "py_dof_manager.hh"
#include "py_dumpable.hh"
#include "py_fe_engine.hh"
#include "py_group_manager.hh"
#include "py_integration_scheme.hh"
#include "py_mesh.hh"
#include "py_model.hh"
#include "py_parser.hh"
#include "py_solver.hh"
#if defined(AKANTU_SOLID_MECHANICS)
#include "py_material.hh"
#include "py_solid_mechanics_model.hh"
#endif
#if defined(AKANTU_DIFFUSION)
#include "py_heat_transfer_model.hh"
#endif
#if defined(AKANTU_COHESIVE_ELEMENT)
#include "py_fragment_manager.hh"
#include "py_solid_mechanics_model_cohesive.hh"
#endif
#if defined(AKANTU_CONTACT_MECHANICS)
#include "py_contact_mechanics_model.hh"
#include "py_model_couplers.hh"
#endif
#if defined(AKANTU_PHASE_FIELD)
#include "py_phase_field_model.hh"
#endif
#if defined(AKANTU_STRUCTURAL_MECHANICS)
#include "py_structural_mechanics_model.hh"
#endif
/* -------------------------------------------------------------------------- */
#include
/* -------------------------------------------------------------------------- */
#include
/* -------------------------------------------------------------------------- */
#include
/* -------------------------------------------------------------------------- */
namespace py = pybind11;
namespace akantu {
void register_all(pybind11::module & mod) {
register_initialize(mod);
register_enums(mod);
register_error(mod);
register_functions(mod);
register_parser(mod);
register_solvers(mod);
- register_group_manager(mod);
register_dumpable(mod);
+ register_group_manager(mod);
register_mesh(mod);
register_fe_engine(mod);
register_integration_schemes(mod);
register_dof_manager(mod);
register_boundary_conditions(mod);
register_model(mod);
register_constitutive_law_selector(mod);
register_constitutive_law_internal_handler(mod);
#if defined(AKANTU_DIFFUSION)
register_heat_transfer_model(mod);
#endif
#if defined(AKANTU_SOLID_MECHANICS)
register_solid_mechanics_model(mod);
register_material(mod);
#endif
#if defined(AKANTU_COHESIVE_ELEMENT)
register_solid_mechanics_model_cohesive(mod);
register_fragment_manager(mod);
#endif
#if defined(AKANTU_STRUCTURAL_MECHANICS)
register_structural_mechanics_model(mod);
#endif
#if defined(AKANTU_CONTACT_MECHANICS)
register_contact_mechanics_model(mod);
register_model_couplers(mod);
#endif
#if defined(AKANTU_PHASE_FIELD)
register_phase_field_model(mod);
register_phase_field_coupler(mod);
#endif
}
} // namespace akantu
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
PYBIND11_MODULE(py11_akantu, mod) {
mod.doc() = "Akantu python interface";
static py::exception akantu_exception(mod,
"Exception");
static py::exception
akantu_exception_nls_not_converged(mod, "NLSNotConvergedException");
py::register_exception_translator([](std::exception_ptr ptr) {
try {
if (ptr) {
std::rethrow_exception(ptr);
}
} catch (akantu::debug::NLSNotConvergedException & e) {
akantu_exception_nls_not_converged(e.info().c_str());
} catch (akantu::debug::Exception & e) {
if (akantu::debug::debugger.printBacktrace()) {
akantu::debug::printBacktrace();
}
akantu_exception(e.info().c_str());
}
});
akantu::register_all(mod);
mod.def("has_mpi",
[]() {
#if defined(AKANTU_USE_MPI)
return true;
#else
return false;
#endif
})
.def("getVersion", &akantu::getVersion);
} // Module akantu
diff --git a/python/py_fe_engine.cc b/python/py_fe_engine.cc
index 4de1f9f0c..71cbcaeab 100644
--- a/python/py_fe_engine.cc
+++ b/python/py_fe_engine.cc
@@ -1,139 +1,211 @@
/**
* Copyright (©) 2019-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* This file is part of Akantu
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see .
*/
/* -------------------------------------------------------------------------- */
#include "py_aka_array.hh"
#include "py_aka_common.hh"
/* -------------------------------------------------------------------------- */
#include
#include
#include
/* -------------------------------------------------------------------------- */
#include
#include
#include
/* -------------------------------------------------------------------------- */
#include
/* -------------------------------------------------------------------------- */
namespace py = pybind11;
/* -------------------------------------------------------------------------- */
namespace akantu {
void register_fe_engine(py::module & mod) {
py::class_(mod, "Element")
.def(py::init([](ElementType type, Int id, GhostType ghost_type) {
return std::make_unique(Element{type, id, ghost_type});
}),
py::arg("type"), py::arg("ghost_type"),
py::arg("ghost_type") = _not_ghost)
.def("__lt__",
[](Element & self, const Element & other) { return (self < other); })
.def("__repr__", [](Element & self) { return std::to_string(self); });
mod.attr("ElementNull") = ElementNull;
py::class_(mod, "FEEngine")
.def(
"getNbIntegrationPoints",
[](FEEngine & fem, ElementType type, GhostType ghost_type) {
return fem.getNbIntegrationPoints(type, ghost_type);
},
py::arg("type"), py::arg("ghost_type") = _not_ghost)
.def("initShapeFunctions", &FEEngine::initShapeFunctions,
py::arg("ghost_type") = _not_ghost)
+
+ .def(
+ "integrate",
+ [](FEEngine & self, const Array & f, Array & intf,
+ Int nb_degree_of_freedom, ElementType type,
+ GhostType ghost_type = _not_ghost,
+ const Array & filter_elements = empty_filter) {
+ self.integrate(f, intf, nb_degree_of_freedom, type, ghost_type,
+ filter_elements);
+ },
+ py::arg("f"), py::arg("intf"), py::arg("nb_degree_of_freedom"),
+ py::arg("type"), py::arg("ghost_type") = _not_ghost,
+ py::arg("filter_elements") = empty_filter)
+ .def(
+ "integrate",
+ [](FEEngine & self, const Array & f, ElementType type,
+ GhostType ghost_type = _not_ghost,
+ const Array & filter_elements = empty_filter) -> Real {
+ return self.integrate(f, type, ghost_type, filter_elements);
+ },
+ py::arg("f"), py::arg("type"), py::arg("ghost_type") = _not_ghost,
+ py::arg("filter_elements") = empty_filter)
.def(
"gradientOnIntegrationPoints",
[](FEEngine & fem, const Array & u, Array & nablauq,
const Int nb_degree_of_freedom, ElementType type,
GhostType ghost_type, const Array & filter_elements) {
fem.gradientOnIntegrationPoints(u, nablauq, nb_degree_of_freedom,
type, ghost_type, filter_elements);
},
py::arg("u"), py::arg("nablauq"), py::arg("nb_degree_of_freedom"),
py::arg("type"), py::arg("ghost_type") = _not_ghost,
py::arg("filter_elements") = empty_filter)
.def(
"interpolateOnIntegrationPoints",
[](FEEngine & self, const Array & u, Array & uq,
Int nb_degree_of_freedom, ElementType type, GhostType ghost_type,
const Array & filter_elements) {
self.interpolateOnIntegrationPoints(
u, uq, nb_degree_of_freedom, type, ghost_type, filter_elements);
},
py::arg("u"), py::arg("uq"), py::arg("nb_degree_of_freedom"),
py::arg("type"), py::arg("ghost_type") = _not_ghost,
py::arg("filter_elements") = empty_filter)
.def(
"interpolateOnIntegrationPoints",
[](FEEngine & self, const Array & u,
std::shared_ptr> uq,
std::shared_ptr> filter_elements) {
self.interpolateOnIntegrationPoints(u, *uq, filter_elements.get());
},
py::arg("u"), py::arg("uq"), py::arg("filter_elements") = nullptr)
+
+ .def(
+ "computeBtD",
+ [](FEEngine & self, const Array & Ds, Array & BtDs,
+ ElementType type, GhostType ghost_type = _not_ghost,
+ const Array & filter_elements = empty_filter) {
+ self.computeBtD(Ds, BtDs, type, ghost_type, filter_elements);
+ },
+ py::arg("Ds"), py::arg("BtDs"), py::arg("type"),
+ py::arg("ghost_type") = _not_ghost,
+ py::arg("filter_elements") = empty_filter)
+
+ .def(
+ "computeBtDB",
+ [](FEEngine & self, const Array & Ds, Array & BtDBs,
+ Int order_d, ElementType type, GhostType ghost_type = _not_ghost,
+ const Array & filter_elements = empty_filter) {
+ self.computeBtDB(Ds, BtDBs, order_d, type, ghost_type,
+ filter_elements);
+ },
+ py::arg("Ds"), py::arg("BtDBs"), py::arg("order_d"), py::arg("type"),
+ py::arg("ghost_type") = _not_ghost,
+ py::arg("filter_elements") = empty_filter)
+
+ .def(
+ "computeNtb",
+ [](FEEngine & self, const Array & bs, Array & Ntbs,
+ ElementType type, GhostType ghost_type = _not_ghost,
+ const Array & filter_elements = empty_filter) {
+ self.computeNtb(bs, Ntbs, type, ghost_type, filter_elements);
+ },
+ py::arg("bs"), py::arg("Ntbs"), py::arg("type"),
+ py::arg("ghost_type") = _not_ghost,
+ py::arg("filter_elements") = empty_filter)
+ .def(
+ "computeNtbN",
+ [](FEEngine & self, const Array & bs, Array & NtbNs,
+ ElementType type, GhostType ghost_type = _not_ghost,
+ const Array & filter_elements = empty_filter) {
+ self.computeNtbN(bs, NtbNs, type, ghost_type, filter_elements);
+ },
+ py::arg("bs"), py::arg("NtbNs"), py::arg("type"),
+ py::arg("ghost_type") = _not_ghost,
+ py::arg("filter_elements") = empty_filter)
.def(
"computeIntegrationPointsCoordinates",
[](FEEngine & self,
std::shared_ptr> coordinates,
std::shared_ptr> filter_elements)
-> decltype(auto) {
return self.computeIntegrationPointsCoordinates(
*coordinates, filter_elements.get());
},
py::arg("coordinates"), py::arg("filter_elements") = nullptr)
.def(
"assembleFieldLumped",
[](FEEngine & fem,
const std::function &, const Element &)> &
field_funct,
const ID & matrix_id, const ID & dof_id, DOFManager & dof_manager,
ElementType type, GhostType ghost_type) {
fem.assembleFieldLumped(field_funct, matrix_id, dof_id, dof_manager,
type, ghost_type);
},
py::arg("field_funct"), py::arg("matrix_id"), py::arg("dof_id"),
py::arg("dof_manager"), py::arg("type"),
py::arg("ghost_type") = _not_ghost)
.def(
"assembleFieldMatrix",
[](FEEngine & fem,
const std::function &, const Element &)> &
field_funct,
const ID & matrix_id, const ID & dof_id, DOFManager & dof_manager,
ElementType type, GhostType ghost_type = _not_ghost) {
fem.assembleFieldMatrix(field_funct, matrix_id, dof_id, dof_manager,
type, ghost_type);
},
py::arg("field_funct"), py::arg("matrix_id"), py::arg("dof_id"),
py::arg("dof_manager"), py::arg("type"),
py::arg("ghost_type") = _not_ghost)
.def("getElementInradius",
[](FEEngine & self, const Element & element) {
return self.getElementInradius(element);
})
.def("getNormalsOnIntegrationPoints",
&FEEngine::getNormalsOnIntegrationPoints, py::arg("type"),
py::arg("ghost_type") = _not_ghost,
- py::return_value_policy::reference);
+ py::return_value_policy::reference)
+ .def("getShapes", &FEEngine::getShapes, py::arg("type"),
+ py::arg("ghost_type") = _not_ghost, py::arg("id") = 0,
+ py::return_value_policy::reference)
+ .def("getShapesDerivatives", &FEEngine::getShapesDerivatives,
+ py::arg("type"), py::arg("ghost_type") = _not_ghost,
+ py::arg("id") = 0, py::return_value_policy::reference);
py::class_(mod, "IntegrationPoint");
}
} // namespace akantu
diff --git a/python/py_group_manager.cc b/python/py_group_manager.cc
index a4d05bb73..5e43fef7a 100644
--- a/python/py_group_manager.cc
+++ b/python/py_group_manager.cc
@@ -1,172 +1,176 @@
/**
* Copyright (©) 2019-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* This file is part of Akantu
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see .
*/
/* -------------------------------------------------------------------------- */
#include "py_aka_array.hh"
/* -------------------------------------------------------------------------- */
+#include
#include
+#include
#include
/* -------------------------------------------------------------------------- */
+#include
+/* -------------------------------------------------------------------------- */
#include
#include
/* -------------------------------------------------------------------------- */
namespace py = pybind11;
/* -------------------------------------------------------------------------- */
namespace akantu {
/* -------------------------------------------------------------------------- */
void register_group_manager(py::module & mod) {
/* ------------------------------------------------------------------------ */
py::class_(mod, "NodeGroup")
.def(
"getNodes",
[](NodeGroup & self) -> decltype(auto) { return self.getNodes(); },
py::return_value_policy::reference)
.def("__len__", &NodeGroup::size)
.def(
"__iter__",
[](const NodeGroup & self) {
return py::make_iterator(self.begin(), self.end());
},
py::keep_alive<0, 1>())
.def("__contains__", [](const NodeGroup & self,
UInt node) { return self.find(node) != -1; })
.def("getName", &NodeGroup::getName)
.def("clear", &NodeGroup::clear)
.def("empty", &NodeGroup::empty)
.def("append", &NodeGroup::append)
.def(
"add",
[](NodeGroup & self, UInt node, bool check_for_duplicate) {
auto && it = self.add(node, check_for_duplicate);
return *it;
},
py::arg("node"), py::arg("check_for_duplicate") = true)
.def("remove", &NodeGroup::add);
/* ------------------------------------------------------------------------ */
- py::class_(mod, "ElementGroup")
+ py::class_(mod, "ElementGroup")
.def(
"getNodeGroup",
[](ElementGroup & self) -> decltype(auto) {
return self.getNodeGroup();
},
py::return_value_policy::reference)
.def("getName", &ElementGroup::getName)
.def(
"getElements",
[](ElementGroup & self) -> decltype(auto) {
return self.getElements();
},
py::return_value_policy::reference)
.def(
"getNodeGroup",
[](ElementGroup & self) -> decltype(auto) {
return self.getNodeGroup();
},
py::return_value_policy::reference)
.def("__len__", [](const ElementGroup & self) { return self.size(); })
.def("clear", [](ElementGroup & self) { self.clear(); })
.def("empty", &ElementGroup::empty)
.def("append", &ElementGroup::append)
.def("optimize", &ElementGroup::optimize)
.def(
"add",
[](ElementGroup & self, const Element & element, bool add_nodes,
bool check_for_duplicate) {
self.add(element, add_nodes, check_for_duplicate);
},
py::arg("element"), py::arg("add_nodes") = false,
py::arg("check_for_duplicate") = true)
.def("fillFromNodeGroup", &ElementGroup::fillFromNodeGroup)
.def("addDimension", &ElementGroup::addDimension);
/* ------------------------------------------------------------------------ */
py::class_(mod, "GroupManager")
.def(
"getElementGroup",
[](GroupManager & self, const std::string & name) -> decltype(auto) {
return self.getElementGroup(name);
},
py::return_value_policy::reference)
.def("iterateElementGroups",
[](GroupManager & self) -> decltype(auto) {
std::vector> groups;
for (auto & group : self.iterateElementGroups()) {
groups.emplace_back(group);
}
return groups;
})
.def("iterateNodeGroups",
[](GroupManager & self) -> decltype(auto) {
std::vector> groups;
for (auto & group : self.iterateNodeGroups()) {
groups.emplace_back(group);
}
return groups;
})
.def("createNodeGroup", &GroupManager::createNodeGroup,
py::return_value_policy::reference)
.def(
"createElementGroup",
[](GroupManager & self, const std::string & id, Int spatial_dimension,
bool replace_group) -> decltype(auto) {
return self.createElementGroup(id, spatial_dimension,
replace_group);
},
py::arg("id"), py::arg("spatial_dimension"),
py::arg("replace_group") = false, py::return_value_policy::reference)
.def("createGroupsFromMeshDataUInt",
&GroupManager::createGroupsFromMeshData)
.def("createElementGroupFromNodeGroup",
&GroupManager::createElementGroupFromNodeGroup, py::arg("name"),
py::arg("node_group"), py::arg("dimension") = _all_dimensions)
.def(
"getNodeGroup",
[](GroupManager & self, const std::string & name) -> decltype(auto) {
return self.getNodeGroup(name);
},
py::return_value_policy::reference)
.def(
"nodeGroups",
[](GroupManager & self) {
std::vector groups;
for (auto & g : self.iterateNodeGroups()) {
groups.push_back(&g);
}
return groups;
},
py::return_value_policy::reference)
.def(
"elementGroups",
[](GroupManager & self) {
std::vector groups;
for (auto & g : self.iterateElementGroups()) {
groups.push_back(&g);
}
return groups;
},
py::return_value_policy::reference)
.def("createBoundaryGroupFromGeometry",
&GroupManager::createBoundaryGroupFromGeometry);
}
} // namespace akantu
diff --git a/python/py_model.cc b/python/py_model.cc
index 21c9a7544..8ea74cd06 100644
--- a/python/py_model.cc
+++ b/python/py_model.cc
@@ -1,202 +1,238 @@
/**
* Copyright (©) 2019-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* This file is part of Akantu
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see .
*/
/* -------------------------------------------------------------------------- */
#include "py_aka_array.hh"
/* -------------------------------------------------------------------------- */
#include
#include
#include
#include
#include
#include
/* -------------------------------------------------------------------------- */
#include
#include
#include
/* -------------------------------------------------------------------------- */
namespace py = pybind11;
/* -------------------------------------------------------------------------- */
namespace akantu {
template
void register_constitutive_law_handler(py::module & mod,
const std::string & name) {
py::class_,
ModelType>(mod, name.c_str(), py::multiple_inheritance())
.def(
"getConstitutiveLaw",
[](ConstitutiveLawsHandler & self,
UInt constitutive_law_id) -> decltype(auto) {
return self.getConstitutiveLaw(constitutive_law_id);
},
py::arg("constitutive_law_id"), py::return_value_policy::reference)
.def(
"getConstitutiveLaw",
[](ConstitutiveLawsHandler & self,
const ID & constitutive_law_name) -> decltype(auto) {
return self.getConstitutiveLaw(constitutive_law_name);
},
py::arg("constitutive_law_name"), py::return_value_policy::reference)
.def(
"getConstitutiveLaw",
[](ConstitutiveLawsHandler & self,
const Element & element) -> decltype(auto) {
return self.getConstitutiveLaw(element);
},
py::arg("element"), py::return_value_policy::reference)
.def("getNbConstitutiveLaws",
&ConstitutiveLawsHandler::getNbConstitutiveLaws)
.def("getConstitutiveLawIndex",
&ConstitutiveLawsHandler::getConstitutiveLawIndex)
.def("setConstitutiveLawSelector",
[](ConstitutiveLawsHandler & self,
std::shared_ptr
constitutive_law_selector) {
self.setConstitutiveLawSelector(constitutive_law_selector);
})
.def("getConstitutiveLawSelector",
&ConstitutiveLawsHandler::getConstitutiveLawSelector)
.def(
"getConstitutiveLawByElement",
[](const ConstitutiveLawsHandler &
self) -> decltype(auto) {
return self.getConstitutiveLawByElement();
},
py::return_value_policy::reference, py::keep_alive<0, 1>())
.def("reassignConstitutiveLaw",
&ConstitutiveLawsHandler::reassignConstitutiveLaw)
.def(
"registerNewConstitutiveLaw",
[](ConstitutiveLawsHandler & self,
const ID & cl_name, const ID & cl_type,
const ID & opt_param) -> decltype(auto) {
return self.registerNewConstitutiveLaw(cl_name, cl_type, opt_param);
},
py::arg("constitutive_law_name"), py::arg("constitutive_law_type"),
py::arg("option") = "", py::return_value_policy::reference)
.def("initConstitutiveLaws",
&ConstitutiveLawsHandler::initConstitutiveLaws)
.def("flattenInternal",
&ConstitutiveLawsHandler::flattenInternal,
py::return_value_policy::reference)
.def("inflateInternal",
&ConstitutiveLawsHandler::inflateInternal,
py::return_value_policy::reference);
}
/* -------------------------------------------------------------------------- */
void register_model(py::module & mod) {
py::class_(mod, "ModelSolver",
py::multiple_inheritance())
.def(
"getNonLinearSolver",
[](ModelSolver & self, const ID & solver_id) -> NonLinearSolver & {
return self.getNonLinearSolver(solver_id);
},
py::arg("solver_id") = "", py::return_value_policy::reference)
.def(
"getTimeStepSolver",
[](ModelSolver & self, const ID & solver_id) -> TimeStepSolver & {
return self.getTimeStepSolver(solver_id);
},
py::arg("solver_id") = "", py::return_value_policy::reference)
.def(
"solveStep",
[](ModelSolver & self, const ID & solver_id) {
self.solveStep(solver_id);
},
py::arg("solver_id") = "")
.def(
"solveStep",
[](ModelSolver & self, SolverCallback & callback,
const ID & solver_id) { self.solveStep(callback, solver_id); },
py::arg("callback"), py::arg("solver_id") = "");
py::class_(mod, "Model", py::multiple_inheritance())
.def("setBaseName", &Model::setBaseName)
.def("setDirectory", &Model::setDirectory)
.def("getFEEngine", &Model::getFEEngine, py::arg("name") = "",
py::return_value_policy::reference)
.def("getFEEngineBoundary", &Model::getFEEngine, py::arg("name") = "",
py::return_value_policy::reference)
.def("addDumpFieldVector", &Model::addDumpFieldVector)
.def("addDumpField", &Model::addDumpField)
.def("setBaseNameToDumper", &Model::setBaseNameToDumper)
.def("addDumpFieldVectorToDumper", &Model::addDumpFieldVectorToDumper)
.def("addDumpFieldToDumper", &Model::addDumpFieldToDumper)
.def("dump", [](Model & self) { self.dump(); })
.def(
"dump", [](Model & self, UInt step) { self.dump(step); },
py::arg("step"))
.def(
"dump",
[](Model & self, Real time, UInt step) { self.dump(time, step); },
py::arg("time"), py::arg("step"))
.def(
"dump",
[](Model & self, const std::string & dumper) { self.dump(dumper); },
py::arg("dumper_name"))
.def(
"dump",
[](Model & self, const std::string & dumper, UInt step) {
self.dump(dumper, step);
},
py::arg("dumper_name"), py::arg("step"))
.def(
"dump",
[](Model & self, const std::string & dumper, Real time, UInt step) {
self.dump(dumper, time, step);
},
py::arg("dumper_name"), py::arg("time"), py::arg("step"))
+
+ .def(
+ "dumpGroup",
+ [](Model & self, const std::string & group_name) {
+ self.dumpGroup(group_name);
+ },
+ py::arg("group_name"))
+ .def(
+ "setGroupDirectory",
+ [](Model & self, const std::string & directory,
+ const std::string & group_name) {
+ self.setGroupDirectory(directory, group_name);
+ },
+ py::arg("directory"), py::arg("group_name"))
+ .def(
+ "setGroupBaseName",
+ [](Model & self, const std::string & basename,
+ const std::string & group_name) {
+ self.setGroupBaseName(basename, group_name);
+ },
+ py::arg("basename"), py::arg("group_name"))
+ .def(
+ "addDumpGroupField",
+ [](Model & self, const std::string & field_id,
+ const std::string & group_name) {
+ self.addDumpGroupField(field_id, group_name);
+ },
+ py::arg("field_id"), py::arg("group_name"))
+
+ .def(
+ "addDumpGroupFieldVector",
+ [](Model & self, const std::string & field_id,
+ const std::string & group_name) {
+ self.addDumpGroupFieldVector(field_id, group_name);
+ },
+ py::arg("field_id"), py::arg("group_name"))
.def("initNewSolver", &Model::initNewSolver)
.def(
"getNewSolver",
[](Model & self, const std::string id,
const TimeStepSolverType & time,
const NonLinearSolverType & type) {
self.getNewSolver(id, time, type);
},
py::return_value_policy::reference)
.def(
"setIntegrationScheme",
[](Model & self, const std::string id, const std::string primal,
const IntegrationSchemeType & scheme_type,
IntegrationScheme::SolutionType solution_type) {
self.setIntegrationScheme(id, primal, scheme_type, solution_type);
},
py::arg("id"), py::arg("primal"), py::arg("scheme_type"),
py::arg("solution_type") =
IntegrationScheme::SolutionType::_not_defined)
.def("getDOFManager", &Model::getDOFManager,
py::return_value_policy::reference)
.def("assembleMatrix", &Model::assembleMatrix);
}
} // namespace akantu
diff --git a/src/fe_engine/fe_engine.hh b/src/fe_engine/fe_engine.hh
index 24bc43b8c..ee6fa7f90 100644
--- a/src/fe_engine/fe_engine.hh
+++ b/src/fe_engine/fe_engine.hh
@@ -1,370 +1,363 @@
/**
* Copyright (©) 2010-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* This file is part of Akantu
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see .
*/
/* -------------------------------------------------------------------------- */
#include "element_type_map.hh"
#include "mesh_events.hh"
/* -------------------------------------------------------------------------- */
#include
/* -------------------------------------------------------------------------- */
#ifndef AKANTU_FE_ENGINE_HH_
#define AKANTU_FE_ENGINE_HH_
namespace akantu {
class Mesh;
class Integrator;
class ShapeFunctions;
class DOFManager;
class Element;
} // namespace akantu
/* -------------------------------------------------------------------------- */
namespace akantu {
/* -------------------------------------------------------------------------- */
/**
* The generic FEEngine class derived in a FEEngineTemplate class
* containing the
* shape functions and the integration method
*/
class FEEngine : public MeshEventHandler {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
FEEngine(Mesh & mesh, Int element_dimension = _all_dimensions,
const ID & id = "fem");
~FEEngine() override;
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
/// pre-compute all the shape functions, their derivatives and the jacobians
virtual void initShapeFunctions(GhostType ghost_type = _not_ghost) = 0;
/// extract the nodal values and store them per element
template
static void
extractNodalToElementField(const Mesh & mesh, const Array & nodal_f,
Array & elemental_f, ElementType type,
GhostType ghost_type = _not_ghost,
const Array & filter_elements = empty_filter);
/// filter a field
template
static void
filterElementalData(const Mesh & mesh, const Array & quad_f,
Array & filtered_f, ElementType type,
GhostType ghost_type = _not_ghost,
const Array & filter_elements = empty_filter);
/* ------------------------------------------------------------------------ */
/* Integration method bridges */
/* ------------------------------------------------------------------------ */
/// integrate f for all elements of type "type"
virtual void
integrate(const Array & f, Array & intf, Int nb_degree_of_freedom,
ElementType type, GhostType ghost_type = _not_ghost,
const Array & filter_elements = empty_filter) const = 0;
/// integrate a scalar value f on all elements of type "type"
[[nodiscard]] virtual Real
integrate(const Array & f, ElementType type,
GhostType ghost_type = _not_ghost,
const Array & filter_elements = empty_filter) const = 0;
- /// integrate f for all integration points of type "type" but don't sum over
- /// all integration points
- virtual void integrateOnIntegrationPoints(
- const Array & f, Array & intf, Int nb_degree_of_freedom,
- ElementType type, GhostType ghost_type = _not_ghost,
- const Array & filter_elements = empty_filter) const = 0;
-
/// integrate one element scalar value on all elements of type "type"
[[nodiscard]] Real integrate(const Ref f,
const Element & element) const {
return integrate(f, element.type, element.element, element.ghost_type);
}
private:
[[nodiscard]] virtual Real
integrate(const Ref f, ElementType type, Idx index,
GhostType ghost_type = _not_ghost) const = 0;
/* ------------------------------------------------------------------------ */
/* compatibility with old FEEngine fashion */
/* ------------------------------------------------------------------------ */
public:
/// get the number of integration points
[[nodiscard]] virtual Int
getNbIntegrationPoints(ElementType type,
GhostType ghost_type = _not_ghost) const = 0;
/// get the precomputed shapes
[[nodiscard]] const virtual Array &
getShapes(ElementType type, GhostType ghost_type = _not_ghost,
Idx id = 0) const = 0;
/// get the derivatives of shapes
[[nodiscard]] virtual const Array &
getShapesDerivatives(ElementType type, GhostType ghost_type = _not_ghost,
Idx id = 0) const = 0;
/// get integration points
[[nodiscard]] virtual const MatrixXr &
getIntegrationPoints(ElementType type,
GhostType ghost_type = _not_ghost) const = 0;
/* ------------------------------------------------------------------------ */
/* Shape method bridges */
/* ------------------------------------------------------------------------ */
/// Compute the gradient nablauq on the integration points of an element type
/// from nodal values u
virtual void gradientOnIntegrationPoints(
const Array & u, Array & nablauq, Int nb_degree_of_freedom,
ElementType type, GhostType ghost_type = _not_ghost,
const Array & filter_elements = empty_filter) const = 0;
/// Interpolate a nodal field u at the integration points of an element type
/// -> uq
virtual void interpolateOnIntegrationPoints(
const Array & u, Array & uq, Int nb_degree_of_freedom,
ElementType type, GhostType ghost_type = _not_ghost,
const Array & filter_elements = empty_filter) const = 0;
/// Interpolate a nodal field u at the integration points of many element
/// types -> uq
virtual void interpolateOnIntegrationPoints(
const Array & u, ElementTypeMapArray & uq,
const ElementTypeMapArray * filter_elements = nullptr) const = 0;
/// pre multiplies a tensor by the shapes derivaties
virtual void
computeBtD(const Array & Ds, Array & BtDs, ElementType type,
GhostType ghost_type = _not_ghost,
const Array & filter_elements = empty_filter) const = 0;
/// left and right multiplies a tensor by the shapes derivaties
virtual void
computeBtDB(const Array & Ds, Array & BtDBs, Int order_d,
ElementType type, GhostType ghost_type = _not_ghost,
const Array & filter_elements = empty_filter) const = 0;
/// left multiples a vector by the shape functions
virtual void
computeNtb(const Array & bs, Array & Ntbs, ElementType type,
GhostType ghost_type = _not_ghost,
const Array & filter_elements = empty_filter) const = 0;
/// left and right multiplies a tensor by the shapes
virtual void
computeNtbN(const Array & bs, Array & NtbNs, ElementType type,
GhostType ghost_type = _not_ghost,
const Array & filter_elements = empty_filter) const = 0;
/// Compute the interpolation point position in the global coordinates for
/// many element types
virtual void computeIntegrationPointsCoordinates(
ElementTypeMapArray & integration_points_coordinates,
const ElementTypeMapArray * filter_elements = nullptr) const = 0;
/// Compute the interpolation point position in the global coordinates for an
/// element type
virtual void computeIntegrationPointsCoordinates(
Array & integration_points_coordinates, ElementType type,
GhostType ghost_type = _not_ghost,
const Array & filter_elements = empty_filter) const = 0;
/// Build pre-computed matrices for interpolation of field form integration
/// points at other given positions (interpolation_points)
virtual void initElementalFieldInterpolationFromIntegrationPoints(
const ElementTypeMapArray & interpolation_points_coordinates,
ElementTypeMapArray & interpolation_points_coordinates_matrices,
ElementTypeMapArray & integration_points_coordinates_inv_matrices,
const ElementTypeMapArray * element_filter) const = 0;
/// interpolate field at given position (interpolation_points) from given
/// values of this field at integration points (field)
virtual void interpolateElementalFieldFromIntegrationPoints(
const ElementTypeMapArray & field,
const ElementTypeMapArray & interpolation_points_coordinates,
ElementTypeMapArray & result, const GhostType ghost_type,
const ElementTypeMapArray * element_filter) const = 0;
/// Interpolate field at given position from given values of this field at
/// integration points (field)
/// using matrices precomputed with
/// initElementalFieldInterplationFromIntegrationPoints
virtual void interpolateElementalFieldFromIntegrationPoints(
const ElementTypeMapArray & field,
const ElementTypeMapArray &
interpolation_points_coordinates_matrices,
const ElementTypeMapArray &
integration_points_coordinates_inv_matrices,
ElementTypeMapArray & result, const GhostType ghost_type,
const ElementTypeMapArray * element_filter) const = 0;
/// interpolate on a phyiscal point inside an element
virtual void interpolate(const Ref real_coords,
const Ref nodal_values,
Ref interpolated,
const Element & element) const = 0;
/// compute the shape on a provided point
virtual void computeShapes(const Ref real_coords, Int elem,
ElementType type, Ref shapes,
GhostType ghost_type = _not_ghost) const = 0;
/// compute the shape derivatives on a provided point
virtual void
computeShapeDerivatives(const Ref real_coords, Int element,
ElementType type, Ref shape_derivatives,
GhostType ghost_type = _not_ghost) const = 0;
/// assembles the lumped version of @f[ \int N^t rho N @f]
virtual void assembleFieldLumped(
const std::function &, const Element &)> & field_funct,
const ID & matrix_id, const ID & dof_id, DOFManager & dof_manager,
ElementType type, GhostType ghost_type = _not_ghost) const = 0;
/// assembles the matrix @f[ \int N^t rho N @f]
virtual void assembleFieldMatrix(
const std::function &, const Element &)> & field_funct,
const ID & matrix_id, const ID & dof_id, DOFManager & dof_manager,
ElementType type, GhostType ghost_type = _not_ghost) const = 0;
/* ------------------------------------------------------------------------ */
/* Other methods */
/* ------------------------------------------------------------------------ */
/// pre-compute normals on integration points
virtual void
computeNormalsOnIntegrationPoints(GhostType ghost_type = _not_ghost) = 0;
/// pre-compute normals on integration points
virtual void
computeNormalsOnIntegrationPoints(const Array & /*field*/,
GhostType /*ghost_type*/ = _not_ghost) {
AKANTU_TO_IMPLEMENT();
}
/// pre-compute normals on integration points
virtual void computeNormalsOnIntegrationPoints(
const Array & /*field*/, Array & /*normal*/,
ElementType /*type*/, GhostType /*ghost_type*/ = _not_ghost) const {
AKANTU_TO_IMPLEMENT();
}
/// function to print the containt of the class
virtual void printself(std::ostream & stream, int indent = 0) const;
private:
/// initialise the class
void init();
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
public:
using ElementTypesIteratorHelper =
ElementTypeMapArray::ElementTypesIteratorHelper;
[[nodiscard]] ElementTypesIteratorHelper
elementTypes(Int dim = _all_dimensions, GhostType ghost_type = _not_ghost,
ElementKind kind = _ek_regular) const;
/// get the dimension of the element handeled by this fe_engine object
AKANTU_GET_MACRO_AUTO(ElementDimension, element_dimension);
/// get the mesh contained in the fem object
AKANTU_GET_MACRO_AUTO(Mesh, mesh);
/// get the mesh contained in the fem object
AKANTU_GET_MACRO_NOT_CONST(Mesh, mesh, Mesh &);
/// get the in-radius of an element
template
[[nodiscard]] static inline Real
getElementInradius(const Eigen::MatrixBase & coord,
ElementType type);
[[nodiscard]] inline Real getElementInradius(const Element & element) const;
/// get the normals on integration points
AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(NormalsOnIntegrationPoints,
normals_on_integration_points, Real);
/// get cohesive element type for a given facet type
[[nodiscard]] static inline constexpr ElementType
getCohesiveElementType(ElementType type_facet);
/// get igfem element type for a given regular type
[[nodiscard]] static inline Vector
getIGFEMElementTypes(ElementType type);
/// get the interpolation element associated to an element type
[[nodiscard]] static inline constexpr auto
getInterpolationType(ElementType el_type);
/// get the shape function class (probably useless: see getShapeFunction in
/// fe_engine_template.hh)
[[nodiscard]] virtual const ShapeFunctions &
getShapeFunctionsInterface() const = 0;
/// get the integrator class (probably useless: see getIntegrator in
/// fe_engine_template.hh)
[[nodiscard]] virtual const Integrator & getIntegratorInterface() const = 0;
AKANTU_GET_MACRO(ID, id, ID);
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
protected:
ID id;
/// spatial dimension of the problem
Int element_dimension;
/// the mesh on which all computation are made
Mesh & mesh;
/// normals at integration points
ElementTypeMapArray normals_on_integration_points;
};
/* -------------------------------------------------------------------------- */
/* inline functions */
/* -------------------------------------------------------------------------- */
/// standard output stream operator
inline std::ostream & operator<<(std::ostream & stream,
const FEEngine & _this) {
_this.printself(stream);
return stream;
}
} // namespace akantu
#include "fe_engine_inline_impl.hh"
#include "fe_engine_template.hh"
#endif /* AKANTU_FE_ENGINE_HH_ */
diff --git a/src/fe_engine/fe_engine_template.hh b/src/fe_engine/fe_engine_template.hh
index b53502f77..a442e0ee7 100644
--- a/src/fe_engine/fe_engine_template.hh
+++ b/src/fe_engine/fe_engine_template.hh
@@ -1,481 +1,474 @@
/**
* Copyright (©) 2010-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* This file is part of Akantu
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see .
*/
/* -------------------------------------------------------------------------- */
#include "fe_engine.hh"
/* -------------------------------------------------------------------------- */
#include
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
#ifndef AKANTU_FE_ENGINE_TEMPLATE_HH_
#define AKANTU_FE_ENGINE_TEMPLATE_HH_
namespace akantu {
class Integrator;
class ShapeFunctions;
} // namespace akantu
namespace akantu {
class DOFManager;
namespace fe_engine {
namespace details {
template struct AssembleLumpedTemplateHelper;
template struct AssembleFieldMatrixHelper;
} // namespace details
} // namespace fe_engine
template struct AssembleFieldMatrixStructHelper;
struct DefaultIntegrationOrderFunctor {
template static inline constexpr int getOrder() {
return ElementClassProperty::polynomial_degree;
}
};
/* -------------------------------------------------------------------------- */
template class I, template class S,
ElementKind kind = _ek_regular,
class IntegrationOrderFunctor = DefaultIntegrationOrderFunctor>
class FEEngineTemplate : public FEEngine {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
using Integ = I;
using Shape = S;
FEEngineTemplate(Mesh & mesh, Int spatial_dimension = _all_dimensions,
const ID & id = "fem", bool do_not_precompute = false);
~FEEngineTemplate() override;
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
/// pre-compute all the shape functions, their derivatives and the jacobians
void initShapeFunctions(GhostType ghost_type = _not_ghost) override;
void initShapeFunctions(const Array & nodes,
GhostType ghost_type = _not_ghost);
/* ------------------------------------------------------------------------ */
/* Integration method bridges */
/* ------------------------------------------------------------------------ */
/// integrate f for all elements of type "type"
void
integrate(const Array & f, Array & intf, Int nb_degree_of_freedom,
ElementType type, GhostType ghost_type = _not_ghost,
const Array & filter_elements = empty_filter) const override;
/// integrate a scalar value on all elements of type "type"
[[nodiscard]] Real
integrate(const Array & f, ElementType type,
GhostType ghost_type = _not_ghost,
const Array & filter_elements = empty_filter) const override;
/// integrate one element scalar value on all elements of type "type"
[[nodiscard]] Real
integrate(const Ref f, ElementType type, Int index,
GhostType ghost_type = _not_ghost) const override;
- /// integrate partially around an integration point (@f$ intf_q = f_q * J_q *
- /// w_q @f$)
- void integrateOnIntegrationPoints(
- const Array & f, Array & intf, Int nb_degree_of_freedom,
- ElementType type, GhostType ghost_type = _not_ghost,
- const Array & filter_elements = empty_filter) const override;
-
private:
template ::value and
kind_ == _ek_regular> * = nullptr>
inline void interpolateImpl(const Eigen::MatrixBase & real_coords,
const Eigen::MatrixBase & nodal_values,
Eigen::MatrixBase & interpolated,
const Element & element) const;
template ::value and
kind_ != _ek_regular> * = nullptr>
inline void interpolateImpl(const Eigen::MatrixBase & /*real_coords*/,
const Eigen::MatrixBase & /*nodal_values*/,
Eigen::MatrixBase & /*interpolated*/,
const Element & /*element*/) const {
AKANTU_TO_IMPLEMENT();
}
public:
/// interpolate on a phyiscal point inside an element
void interpolate(const Ref real_coords,
const Ref nodal_values,
Ref interpolated,
const Element & element) const override;
/// get the number of integration points
[[nodiscard]] Int
getNbIntegrationPoints(ElementType type,
GhostType ghost_type = _not_ghost) const override;
/// get shapes precomputed
[[nodiscard]] const Array & getShapes(ElementType type,
GhostType ghost_type = _not_ghost,
Int id = 0) const override;
/// get the derivatives of shapes
[[nodiscard]] const Array &
getShapesDerivatives(ElementType type, GhostType ghost_type = _not_ghost,
Int id = 0) const override;
/// get integration points
[[nodiscard]] inline const Matrix &
getIntegrationPoints(ElementType type,
GhostType ghost_type = _not_ghost) const override;
/* ------------------------------------------------------------------------ */
/* Shape method bridges */
/* ------------------------------------------------------------------------ */
/// compute the gradient of a nodal field on the integration points
void gradientOnIntegrationPoints(
const Array & u, Array & nablauq, Int nb_degree_of_freedom,
ElementType type, GhostType ghost_type = _not_ghost,
const Array & filter_elements = empty_filter) const override;
/// interpolate a nodal field on the integration points
void interpolateOnIntegrationPoints(
const Array & u, Array & uq, Int nb_degree_of_freedom,
ElementType type, GhostType ghost_type = _not_ghost,
const Array & filter_elements = empty_filter) const override;
/// interpolate a nodal field on the integration points given a
/// by_element_type
void interpolateOnIntegrationPoints(
const Array & u, ElementTypeMapArray & uq,
const ElementTypeMapArray * filter_elements =
nullptr) const override;
/// pre multiplies a tensor by the shapes derivaties
void
computeBtD(const Array & Ds, Array & BtDs, ElementType type,
GhostType ghost_type,
const Array & filter_elements = empty_filter) const override;
/// left and right multiplies a tensor by the shapes derivaties
void
computeBtDB(const Array & Ds, Array & BtDBs, Int order_d,
ElementType type, GhostType ghost_type,
const Array & filter_elements = empty_filter) const override;
/// left multiples a vector by the shape functions
void computeNtb(const Array & bs, Array & Ntbs, ElementType type,
GhostType ghost_type,
const Array & filter_elements) const override;
/// left and right multiplies a tensor by the shapes
void
computeNtbN(const Array & bs, Array & NtbNs, ElementType type,
GhostType ghost_type,
const Array & filter_elements = empty_filter) const override;
/// compute the position of integration points given by an element_type_map
/// from nodes position
inline void computeIntegrationPointsCoordinates(
ElementTypeMapArray & quadrature_points_coordinates,
const ElementTypeMapArray * filter_elements =
nullptr) const override;
/// compute the position of integration points from nodes position
inline void computeIntegrationPointsCoordinates(
Array & quadrature_points_coordinates, ElementType type,
GhostType ghost_type = _not_ghost,
const Array & filter_elements = empty_filter) const override;
/// interpolate field at given position (interpolation_points) from given
/// values of this field at integration points (field)
inline void interpolateElementalFieldFromIntegrationPoints(
const ElementTypeMapArray & field,
const ElementTypeMapArray & interpolation_points_coordinates,
ElementTypeMapArray & result, GhostType ghost_type,
const ElementTypeMapArray * element_filter) const override;
/// Interpolate field at given position from given values of this field at
/// integration points (field)
/// using matrices precomputed with
/// initElementalFieldInterplationFromIntegrationPoints
inline void interpolateElementalFieldFromIntegrationPoints(
const ElementTypeMapArray & field,
const ElementTypeMapArray &
interpolation_points_coordinates_matrices,
const ElementTypeMapArray & quad_points_coordinates_inv_matrices,
ElementTypeMapArray & result, GhostType ghost_type,
const ElementTypeMapArray * element_filter) const override;
/// Build pre-computed matrices for interpolation of field form integration
/// points at other given positions (interpolation_points)
inline void initElementalFieldInterpolationFromIntegrationPoints(
const ElementTypeMapArray & interpolation_points_coordinates,
ElementTypeMapArray & interpolation_points_coordinates_matrices,
ElementTypeMapArray & quad_points_coordinates_inv_matrices,
const ElementTypeMapArray * element_filter = nullptr) const override;
/// find natural coords from real coords provided an element
void inverseMap(const Ref real_coords, Int element,
ElementType type, Ref natural_coords,
GhostType ghost_type = _not_ghost) const;
/// return true if the coordinates provided are inside the element, false
/// otherwise
inline bool contains(const Vector & real_coords, Int element,
ElementType type,
GhostType ghost_type = _not_ghost) const;
private:
template ::value and
kind_ != _ek_cohesive> * = nullptr>
inline void computeShapesImpl(const Eigen::MatrixBase & real_coords,
Int element, ElementType type,
Eigen::MatrixBase & shapes,
GhostType ghost_type = _not_ghost) const;
template ::value and
kind_ == _ek_cohesive> * = nullptr>
inline void computeShapesImpl(const Eigen::MatrixBase & /*real_coords*/,
Int /*element*/, ElementType /*type*/,
Eigen::MatrixBase & /*shapes*/,
GhostType /*ghost_type*/ = _not_ghost) const {
AKANTU_TO_IMPLEMENT();
}
template and kind_ != _ek_cohesive> * =
nullptr>
inline void
computeShapeDerivativesImpl(const Eigen::MatrixBase & real_coords,
Int element, ElementType type,
Eigen::MatrixBase & shape_derivatives,
GhostType ghost_type = _not_ghost) const;
template and kind_ == _ek_cohesive> * =
nullptr>
inline void
computeShapeDerivativesImpl(const Eigen::MatrixBase & /*real_coords*/,
Int /*element*/, ElementType /*type*/,
Eigen::MatrixBase & /*shape_derivatives*/,
GhostType /*ghost_type*/ = _not_ghost) const {
AKANTU_TO_IMPLEMENT();
}
public:
/// compute the shape on a provided point
inline void computeShapes(const Ref real_coords, Int element,
ElementType type, Ref shapes,
GhostType ghost_type = _not_ghost) const override {
this->template computeShapesImpl(real_coords, element, type, shapes,
ghost_type);
}
/// compute the shape derivatives on a provided point
inline void
computeShapeDerivatives(const Ref