diff --git a/CMakeLists.txt b/CMakeLists.txt
index 5ef7f8c6ee0c51cadbe6049b1aad59ae641c1288..e51f09d5a1e1f2a958c9a3093fe1275db6337b12 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -225,6 +225,7 @@ set(PUGS_ENABLE_COSTO OFF CACHE STRING
if (PUGS_ENABLE_COSTO MATCHES "^(AUTO|ON)$")
find_package(Costo 0.0.1 COMPONENTS CPO)
+ find_package(Python3 3.10 EXACT COMPONENTS Interpreter Development NumPy)
if(Costo_FOUND)
else()
@@ -239,6 +240,8 @@ endif()
if (${Costo_FOUND})
include_directories(SYSTEM ${Costo_INCLUDE_DIR})
+ include_directories(SYSTEM ${Python3_INCLUDE_DIRS})
+ include_directories(SYSTEM ${Python3_NumPy_INCLUDE_DIRS})
link_directories(${Costo_LIBRARY_DIR})
else()
if (PUGS_ENABLE_COSTO MATCHES "^ON$")
diff --git a/src/language/modules/CouplageModule.cpp b/src/language/modules/CouplageModule.cpp
index 3f53aa49a7bb067488df319e02abc650f9fffe18..f4cb22ab12919df5faabe262266e23b167ed0ccc 100644
--- a/src/language/modules/CouplageModule.cpp
+++ b/src/language/modules/CouplageModule.cpp
@@ -50,8 +50,18 @@ CouplageModule::CouplageModule()
const std::shared_ptr<const MeshVariant> mesh_v,
const int64_t location,
const int64_t dest, const int64_t tag) -> void {
- return Serializer().BC(boundary, mesh_v, location, dest, tag);
- }));
+ return Serializer().BC(boundary, mesh_v, location, dest, tag);
+ }));
+
+
+ this->_addBuiltinFunction("cpl_build_graph", std::function([](const std::shared_ptr<const MeshVariant> mesh_v,
+ const std::shared_ptr<const IBoundaryDescriptor>& boundary) -> void {
+ return Serializer().build_graph(mesh_v, boundary);
+ }));
+
+ this->_addBuiltinFunction("cpl_get_dt", std::function([](const double& dt) -> double {
+ return Serializer().cpl_get_dt(dt);
+ }));
this->_addBuiltinFunction("serialize_field",
std::function([](const std::shared_ptr<const MeshVariant> mesh,
diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index 3526c4d1d81c35fd7e2dd92b6c86867a09c83670..29932e8ab1b5bad310ec7d1ef8551968ee9b28d8 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -17,7 +17,6 @@
#include <mesh/MeshDataManager.hpp>
#include <mesh/MeshRandomizer.hpp>
#include <mesh/MeshSmoother.hpp>
-// #include <mesh/NumberedBoundaryDescriptor.hpp>
#include <mesh/MeshTraits.hpp>
#include <scheme/AcousticSolver.hpp>
#include <scheme/AxisBoundaryConditionDescriptor.hpp>
diff --git a/src/scheme/HyperelasticSolver.cpp b/src/scheme/HyperelasticSolver.cpp
index ea1b0f525beb3cea6f29e29c460610cda2290258..b8679dbe2185bde0c3b1244b276831aaf86efb52 100644
--- a/src/scheme/HyperelasticSolver.cpp
+++ b/src/scheme/HyperelasticSolver.cpp
@@ -408,7 +408,7 @@ class HyperelasticSolverHandler::HyperelasticSolver final : public HyperelasticS
#endif // PUGS_HAS_COSTO
void _applySymmetryBC(const BoundaryConditionList& bc_list, NodeValue<Rdxd>& Ar, NodeValue<Rd>& br) const;
void _applyVelocityBC(const BoundaryConditionList& bc_list, NodeValue<Rdxd>& Ar, NodeValue<Rd>& br) const;
- void _applyCouplingBC(const BoundaryConditionList& bc_list, NodeValue<Rdxd>& Ar, NodeValue<Rd>& br) const;
+ void _applyCouplingBC(const BoundaryConditionList& bc_list, const MeshType& mesh, NodeValue<Rdxd>& Ar, NodeValue<Rd>& br) const;
void _applyCouplingBC2(const MeshType& mesh,
const BoundaryConditionList& bc_list,
@@ -498,7 +498,7 @@ class HyperelasticSolverHandler::HyperelasticSolver final : public HyperelasticS
const BoundaryConditionList bc_list = this->_getBCList(mesh, bc_descriptor_list);
this->_applyBoundaryConditions(bc_list, mesh, Ar, br);
- this->_applyCouplingBC(bc_list, Ar, br);
+ this->_applyCouplingBC(bc_list, mesh, Ar, br);
synchronize(Ar);
synchronize(br);
@@ -864,7 +864,8 @@ HyperelasticSolverHandler::HyperelasticSolver<MeshType>::_applyVelocityBC(const
template <MeshConcept MeshType>
void
HyperelasticSolverHandler::HyperelasticSolver<MeshType>::_applyCouplingBC(const BoundaryConditionList& bc_list,
- NodeValue<Rdxd>& Ar,
+ const MeshType& mesh,
+ NodeValue<Rdxd>& Ar,
NodeValue<Rd>& br) const
{
for (const auto& boundary_condition : bc_list) {
@@ -874,68 +875,102 @@ HyperelasticSolverHandler::HyperelasticSolver<MeshType>::_applyCouplingBC(const
if constexpr (std::is_same_v<CouplingBoundaryCondition, T>) {
const auto& node_list = bc.nodeList();
- /* std::cout << "\033[01;31m" */
- /* << "un applyCoupling" */
- /* << "Dimension: " << Dimension << "\033[00;00m" << std::endl; */
- /* std::cout << "\033[01;31m" */
- /* << "node_list.size()" << node_list.size() << "\033[00;00m" << std::endl; */
-
#ifdef PUGS_HAS_COSTO
- auto costo = parallel::Messenger::getInstance().myCoupling;
- const int dest = costo->myGlobalSize() - 1;
- int tag = 200;
+ auto costo = parallel::Messenger::getInstance().myPyEmbeding;
+ const auto node_is_owned = mesh.connectivity().nodeIsOwned().arrayView();
+
+
+ // auto costo = parallel::Messenger::getInstance().myCoupling;
+ // const int dest = costo->myGlobalSize() - 1;
+ // int tag = 200;
+
+ // TO DO COMPUTE ONLY ONCE AT THE BEGINING
+ size_t sz = 0;
+ for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
+ const NodeId r = node_list[i_node];
+ if (node_is_owned[r]){
+ sz += 1;
+ }
+ }
+ Array<TinyVector<Dimension>> br_extracted(node_list.size());
+ Array<TinyMatrix<Dimension>> Ar_extracted(node_list.size());
+
+
std::vector<int> shape;
shape.resize(3);
- shape[0] = node_list.size();
+ shape[0] = sz;
shape[1] = Dimension + Dimension * Dimension;
/* shape[1] = Dimension; */
shape[2] = 0;
std::vector<double> data;
data.resize(shape[0] * shape[1] + shape[2]);
- Array<TinyVector<Dimension>> br_extracted(node_list.size());
- Array<TinyMatrix<Dimension>> Ar_extracted(node_list.size());
+
parallel_for(
node_list.size(), PUGS_LAMBDA(size_t i_node) {
- NodeId node_id = node_list[i_node];
- for (size_t i_dim = 0; i_dim < Dimension; i_dim++) {
- br_extracted[i_node] = br[node_id];
- Ar_extracted[i_node] = Ar[node_id];
- }
+ const NodeId node_id = node_list[i_node];
+ if (node_is_owned[node_id]){
+ for (size_t i_dim = 0; i_dim < Dimension; i_dim++) {
+ br_extracted[i_node] = br[node_id];
+ Ar_extracted[i_node] = Ar[node_id];
+ }
+ }
});
+
+ size_t ct_node = 0;
/*TODO: serializer directement Ar et br pour eviter une copie*/
for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
- for (unsigned short i_dim = 0; i_dim < Dimension; ++i_dim) {
- data[i_node * shape[1] + i_dim] = br_extracted[i_node][i_dim];
- for (size_t j_dim = 0; j_dim < Dimension; j_dim++) {
- data[i_node * shape[1] + Dimension + i_dim * Dimension + j_dim] = Ar_extracted[i_node](i_dim, j_dim);
- }
- }
+ const NodeId r = node_list[i_node];
+ if (node_is_owned[r]){
+ for (unsigned short i_dim = 0; i_dim < Dimension; ++i_dim) {
+ data[ct_node * shape[1] + i_dim] = br_extracted[i_node][i_dim];
+ for (size_t j_dim = 0; j_dim < Dimension; j_dim++) {
+ data[ct_node * shape[1] + Dimension + i_dim * Dimension + j_dim] = Ar_extracted[i_node](i_dim, j_dim);
+ }
+ }
+ ct_node += 1;
+ }
}
- costo->sendData(shape, data, dest, tag);
- std::vector<double> recv;
- tag = 210;
- costo->recvData(recv, dest, tag);
+ costo->Contexts->populate_context_old("cpl_bc", "shape", shape);
+ costo->Contexts->populate_context_old("cpl_bc", "Ar_Br", data);
+ costo->call_py_function("cpl_bc");
+ // costo->sendData(shape, data, dest, tag);
+ // std::vector<double> recv;
+
+ // tag = 210;
+ // costo->recvData(recv, dest, tag);
+ ct_node = 0;
for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
- for (unsigned short i_dim = 0; i_dim < Dimension; ++i_dim) {
- br_extracted[i_node][i_dim] = recv[i_node * shape[1] + i_dim];
- for (size_t j_dim = 0; j_dim < Dimension; j_dim++) {
- Ar_extracted[i_node](i_dim, j_dim) = recv[i_node * shape[1] + Dimension + i_dim * Dimension + j_dim];
- }
- }
+ const NodeId r = node_list[i_node];
+ if (node_is_owned[r]){
+ for (unsigned short i_dim = 0; i_dim < Dimension; ++i_dim) {
+ br_extracted[i_node][i_dim] = data[ct_node * shape[1] + i_dim];
+ for (size_t j_dim = 0; j_dim < Dimension; j_dim++) {
+ Ar_extracted[i_node](i_dim, j_dim) = data[ct_node * shape[1] + Dimension + i_dim * Dimension + j_dim];
+ }
+ }
+ ct_node += 1;
+ }
+
+ // std::cout<<"Br: "<<br_extracted[i_node]<<std::endl;
+ // std::cout<<"Ar: "<<Ar_extracted[i_node]<<std::endl;
}
- costo->m = static_cast<int>(recv[shape[0] * shape[1] + shape[2]]);
+ // costo->m = static_cast<int>(recv[shape[0] * shape[1] + shape[2]]);
+
+ // std::cout << "\033[01;32m"
+ // << "q recv: " << costo->m << "\033[00;00m" << std::endl;
+
- std::cout << "\033[01;32m"
- << "q recv: " << costo->m << "\033[00;00m" << std::endl;
parallel_for(
node_list.size(), PUGS_LAMBDA(size_t i_node) {
- NodeId node_id = node_list[i_node];
- for (size_t i_dim = 0; i_dim < Dimension; i_dim++) {
- br[node_id] = br_extracted[i_node];
- Ar[node_id] = Ar_extracted[i_node];
- }
+ const NodeId node_id = node_list[i_node];
+ if (node_is_owned[node_id]){
+ for (size_t i_dim = 0; i_dim < Dimension; i_dim++) {
+ br[node_id] = br_extracted[i_node];
+ Ar[node_id] = Ar_extracted[i_node];
+ }
+ }
});
#endif // PUGS_HAS_COSTO
@@ -1308,6 +1343,8 @@ class HyperelasticSolverHandler::HyperelasticSolver<MeshType>::CouplingBoundaryC
CouplingBoundaryCondition(const MeshNodeBoundary& mesh_node_boundary)
: m_mesh_node_boundary{mesh_node_boundary}
{
+
+ // std::cout<<"Build coupling Boundary condition"<<std::endl;
;
}
diff --git a/src/utils/Messenger.cpp b/src/utils/Messenger.cpp
index e8262c55457f46b0ec0a387afb580eafaed7b448..2c2466961a26d712d29c87dcedb5498560d93f62 100644
--- a/src/utils/Messenger.cpp
+++ b/src/utils/Messenger.cpp
@@ -35,8 +35,10 @@ Messenger::Messenger([[maybe_unused]] int& argc, [[maybe_unused]] char* argv[])
#ifdef PUGS_HAS_COSTO
if (CommunicatorManager::hasSplitColor()) {
- myCoupling = new costo::coupling(CommunicatorManager::splitColor());
+ myPyEmbeding = new costo::embedingPython();
+ myPyEmbeding->call_py_function("cpl_initialize");
+ myCoupling = new costo::coupling(CommunicatorManager::splitColor());
m_pugs_comm_world = myCoupling->getLocalComm();
/* auto res = */
/* MPI_Comm_split(MPI_COMM_WORLD, CommunicatorManager::splitColor(), key, &m_pugs_comm_world); */
diff --git a/src/utils/Messenger.hpp b/src/utils/Messenger.hpp
index d551fda3efbcd2eb11edc2e0583c223fe2b3ff29..9783efffa974383e47445fa80b11226b2871603b 100644
--- a/src/utils/Messenger.hpp
+++ b/src/utils/Messenger.hpp
@@ -18,6 +18,7 @@
#ifdef PUGS_HAS_COSTO
#include <CPO/couplingObjects.h>
+#include <CPO/embedingPython.h>
#endif // PUGS_HAS_COSTO
#include <utils/CommunicatorManager.hpp>
@@ -413,6 +414,8 @@ class Messenger
public:
#ifdef PUGS_HAS_COSTO
costo::coupling* myCoupling = nullptr;
+ costo::embedingPython* myPyEmbeding;
+
#endif // PUGS_HAS_COSTO
static void create(int& argc, char* argv[]);
static void destroy();
diff --git a/src/utils/Serializer.cpp b/src/utils/Serializer.cpp
index 9227e2ec7c3e629afce9e4afb9d56c0549fcfcff..a28b5a20c5dd38e7ec4f782d1a712de040913314 100644
--- a/src/utils/Serializer.cpp
+++ b/src/utils/Serializer.cpp
@@ -11,7 +11,6 @@ Serializer::mesh(const std::shared_ptr<const MeshVariant> mesh_v,
if (not mesh_v) {
throw NormalError("discrete functions are not defined on the same mesh");
}
-
std::visit(
[&](auto&& p_mesh) {
using MeshType = mesh_type_t<decltype(p_mesh)>;
@@ -62,6 +61,82 @@ Serializer::mesh(const std::shared_ptr<const MeshVariant> mesh_v,
mesh_v->variant());
}
+double Serializer::cpl_get_dt(const double& dt)
+{
+ std::vector<double> new_dt;
+ new_dt.resize(1);
+ new_dt[0] = dt*1.0;
+
+ auto costo = parallel::Messenger::getInstance().myPyEmbeding;
+ costo->Contexts->populate_context_old("cpl_get_dt", "dt", new_dt);
+ costo->call_py_function("cpl_get_dt");
+ return new_dt[0];
+}
+
+
+void
+Serializer::build_graph(const std::shared_ptr<const MeshVariant> mesh_v,
+ const IBoundaryDescriptor& boundary)
+{
+
+ if (not mesh_v) {
+ throw NormalError("discrete functions are not defined on the same mesh");
+ }
+ std::visit(
+ [&](auto&& p_mesh) {
+ using MeshType = mesh_type_t<decltype(p_mesh)>;
+ static const auto Dimension = MeshType::Dimension;
+ using Rd = TinyVector<Dimension>;
+ if constexpr (is_polygonal_mesh_v<MeshType>) {
+ auto costo = parallel::Messenger::getInstance().myPyEmbeding;
+ std::vector<int> shape;
+ std::vector<double> pts;
+ const std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(p_mesh);
+
+ const auto xr = mesh->xr();
+
+ const auto& connectivity = mesh->connectivity();
+ const auto node_is_owned = connectivity.nodeIsOwned().arrayView();
+
+ MeshNodeBoundary mesh_node_boundary = getMeshNodeBoundary(*mesh, boundary);
+ const auto node_list = mesh_node_boundary.nodeList();
+
+ Array<TinyVector<Dimension>> positions(node_list.size());
+
+ size_t sz = 0;
+ for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
+ const NodeId r = node_list[i_node];
+ std::cout<<"Serialize node id: "<<r<<std::endl;
+ if (node_is_owned[r]){
+ for (unsigned short i = 0; i < Dimension; ++i) {
+ positions[sz][i] = xr[r][i];
+ }
+ sz += 1;
+ }
+ }
+
+ /*TODO: serializer directement position pour eviter une copie*/
+ pts.resize(sz * Dimension);
+ for (unsigned short r = 0; r < sz; ++r) {
+ for (unsigned short j = 0; j < Dimension; ++j) {
+ pts[Dimension * r + j] = positions[r][j];
+ }
+ }
+ shape.resize(3);
+ shape[0] = sz;
+ shape[1] = Dimension;
+ shape[2] = 0;
+
+ costo->Contexts->populate_context_old("cpl_build_graph", "shape", shape);
+ costo->Contexts->populate_context_old("cpl_build_graph", "mesh", pts);
+ costo->call_py_function("cpl_build_graph");
+ } else {
+ throw NormalError("unexpected mesh type");
+ }
+ },
+ mesh_v->variant());
+}
+
void Serializer::interfaces(const std::shared_ptr<const MeshVariant> mesh_v,
const std::vector<std::shared_ptr<const IInterfaceDescriptor>>& i_interface_list,
const int64_t dest,
diff --git a/src/utils/Serializer.hpp b/src/utils/Serializer.hpp
index 0060afd7937be5862d2bb4bfe4561b39c3644615..20bd6a9980ff4b829ae2e81c2f6b1147858bfab2 100644
--- a/src/utils/Serializer.hpp
+++ b/src/utils/Serializer.hpp
@@ -59,6 +59,20 @@ public:
const int64_t dest,
const int64_t tag);
+
+
+ void build_graph(const std::shared_ptr<const MeshVariant> mesh_v,
+ const std::shared_ptr<const IBoundaryDescriptor>& boundary)
+ {
+ this->build_graph(mesh_v,*boundary);
+ }
+
+ void build_graph(const std::shared_ptr<const MeshVariant> mesh_v,
+ const IBoundaryDescriptor& boundary);
+
+ double cpl_get_dt(const double& dt);
+
+
void
BC(const std::shared_ptr<const IBoundaryDescriptor>& boundary,
const std::shared_ptr<const MeshVariant> mesh_v,