diff --git a/src/mesh/ImplicitMeshSmoother.cpp b/src/mesh/ImplicitMeshSmoother.cpp
index 43538914a2069268254c8d76d765bdf0ce39d7bc..b73d6ec3d425caaf7c0bd327e6f6b554fa74df0d 100644
--- a/src/mesh/ImplicitMeshSmoother.cpp
+++ b/src/mesh/ImplicitMeshSmoother.cpp
@@ -79,46 +79,49 @@ class ImplicitMeshSmootherHandler::ImplicitMeshSmoother
}
void
- _browseBC(NodeValue<bool>& is_fixed) const
+ _browseBC(NodeValue<bool>& is_fixed, NodeValue<int>& is_symmetric) const
{
+ int n_sym = 1;
for (auto&& boundary_condition : m_boundary_condition_list) {
std::visit(
[&](auto&& bc) {
using BCType = std::decay_t<decltype(bc)>;
- // if constexpr (std::is_same_v<BCType, SymmetryBoundaryCondition>) {
- // const Rd& n = bc.outgoingNormal();
+ if constexpr (std::is_same_v<BCType, SymmetryBoundaryCondition>) {
+ // const Rd& n = bc.outgoingNormal();
- // const Rdxd I = identity;
- // const Rdxd nxn = tensorProduct(n, n);
- // const Rdxd P = I - nxn;
+ // const Rdxd I = identity;
+ // const Rdxd nxn = tensorProduct(n, n);
+ // const Rdxd P = I - nxn;
- // const Array<const NodeId>& node_list = bc.nodeList();
- // parallel_for(
- // node_list.size(), PUGS_LAMBDA(const size_t i_node) {
- // const NodeId node_id = node_list[i_node];
-
- // shift[node_id] = P * shift[node_id];
- // });
-
- // } else if constexpr (std::is_same_v<BCType, AxisBoundaryCondition>) {
- // if constexpr (Dimension > 1) {
- // const Rd& t = bc.direction();
+ const Array<const NodeId>& node_list = bc.nodeList();
+ parallel_for(
+ node_list.size(), PUGS_LAMBDA(const size_t i_node) {
+ const NodeId node_id = node_list[i_node];
+ if (is_symmetric[node_id] != 0) {
+ is_fixed[node_id] = true;
+ } else if (!is_fixed[node_id]) {
+ is_symmetric[node_id] = n_sym;
+ }
+ });
+ n_sym += 1;
+ // } else if constexpr (std::is_same_v<BCType, AxisBoundaryCondition>) {
+ // if constexpr (Dimension > 1) {
+ // const Rd& t = bc.direction();
- // const Rdxd txt = tensorProduct(t, t);
+ // const Rdxd txt = tensorProduct(t, t);
- // const Array<const NodeId>& node_list = bc.nodeList();
- // parallel_for(
- // node_list.size(), PUGS_LAMBDA(const size_t i_node) {
- // const NodeId node_id = node_list[i_node];
+ // const Array<const NodeId>& node_list = bc.nodeList();
+ // parallel_for(
+ // node_list.size(), PUGS_LAMBDA(const size_t i_node) {
+ // const NodeId node_id = node_list[i_node];
- // shift[node_id] = txt * shift[node_id];
- // });
- // } else {
- // throw UnexpectedError("AxisBoundaryCondition make no sense in dimension 1");
- // }
+ // shift[node_id] = txt * shift[node_id];
+ // });
+ // } else {
+ // throw UnexpectedError("AxisBoundaryCondition make no sense in dimension 1");
+ // }
- // } else
- if constexpr (std::is_same_v<BCType, FixedBoundaryCondition>) {
+ } else if constexpr (std::is_same_v<BCType, FixedBoundaryCondition>) {
const Array<const NodeId>& node_list = bc.nodeList();
parallel_for(
node_list.size(), PUGS_LAMBDA(const size_t i_node) {
@@ -167,22 +170,23 @@ class ImplicitMeshSmootherHandler::ImplicitMeshSmoother
int local_free_id = 0;
int local_fixed_id = 0;
non_zeros_free.fill(1);
+ // Warning this is not correct
non_zeros_fixed.fill(0);
for (NodeId n_id = 0; n_id < m_given_mesh.numberOfNodes(); n_id++) {
- // std::cout << " n_id " << n_id << "\n";
if (is_fixed[n_id]) {
gid_node_fixed[local_fixed_id] = n_id;
- for (size_t i_face = 0; i_face < node_to_face_matrix[n_id].size(); ++i_face) {
- FaceId face_id = node_to_face_matrix[n_id][i_face];
- for (size_t i_node = 0; i_node < face_to_node_matrix[face_id].size(); ++i_node) {
- NodeId node_id = face_to_node_matrix[face_id][i_node];
- if ((node_id == n_id) or (is_fixed[node_id])) {
- continue;
- } else {
- non_zeros_fixed[local_fixed_id] += 1;
- }
- }
- }
+ // for (size_t i_face = 0; i_face < node_to_face_matrix[n_id].size(); ++i_face) {
+ // FaceId face_id = node_to_face_matrix[n_id][i_face];
+
+ // for (size_t i_node = 0; i_node < face_to_node_matrix[face_id].size(); ++i_node) {
+ // NodeId node_id = face_to_node_matrix[face_id][i_node];
+ // if ((node_id == n_id) or (is_fixed[node_id])) {
+ // continue;
+ // } else {
+ // non_zeros_fixed[local_fixed_id] += 1;
+ // }
+ // }
+ // }
local_fixed_id++;
} else {
gid_node_free[local_free_id] = n_id;
@@ -190,8 +194,10 @@ class ImplicitMeshSmootherHandler::ImplicitMeshSmoother
FaceId face_id = node_to_face_matrix[n_id][i_face];
for (size_t i_node = 0; i_node < face_to_node_matrix[face_id].size(); ++i_node) {
NodeId node_id = face_to_node_matrix[face_id][i_node];
- if ((node_id == n_id) or (is_fixed[node_id])) {
+ if ((node_id == n_id)) {
continue;
+ } else if ((is_fixed[node_id])) {
+ non_zeros_fixed[local_free_id] += 1;
} else {
non_zeros_free[local_free_id] += 1;
}
@@ -205,6 +211,7 @@ class ImplicitMeshSmootherHandler::ImplicitMeshSmoother
_fillMatrix(CRSMatrixDescriptor<double>& Afree,
CRSMatrixDescriptor<double>& Afixed,
const NodeValue<bool>& is_fixed,
+ const NodeValue<int>& is_symmetric,
const NodeValue<int>& id_free,
const NodeValue<int>& id_fixed) const
// ,
@@ -219,18 +226,37 @@ class ImplicitMeshSmootherHandler::ImplicitMeshSmoother
continue;
} else {
Assert(id_free[n_id] == local_free_id, "bad matrix definition");
- for (size_t i_face = 0; i_face < node_to_face_matrix[n_id].size(); ++i_face) {
- FaceId face_id = node_to_face_matrix[n_id][i_face];
- for (size_t i_node = 0; i_node < face_to_node_matrix[face_id].size(); ++i_node) {
- NodeId node_id = face_to_node_matrix[face_id][i_node];
- if (node_id == n_id) {
- continue;
- } else if ((is_fixed[node_id])) {
- Afree(local_free_id, local_free_id) += 1;
- Afixed(local_free_id, id_fixed[node_id]) -= 1;
- } else {
- Afree(local_free_id, local_free_id) += 1;
- Afree(local_free_id, id_free[node_id]) -= 1;
+ if (is_symmetric[n_id]) {
+ for (size_t i_face = 0; i_face < node_to_face_matrix[n_id].size(); ++i_face) {
+ FaceId face_id = node_to_face_matrix[n_id][i_face];
+ for (size_t i_node = 0; i_node < face_to_node_matrix[face_id].size(); ++i_node) {
+ NodeId node_id = face_to_node_matrix[face_id][i_node];
+ if (node_id == n_id) {
+ continue;
+ } else if ((is_fixed[node_id])) {
+ Afree(local_free_id, local_free_id) += 1;
+ Afixed(local_free_id, id_fixed[node_id]) -= 1;
+ } else if (is_symmetric[node_id] == is_symmetric[n_id]) {
+ Afree(local_free_id, local_free_id) += 1;
+ Afree(local_free_id, id_free[node_id]) -= 1;
+ }
+ }
+ }
+
+ } else {
+ for (size_t i_face = 0; i_face < node_to_face_matrix[n_id].size(); ++i_face) {
+ FaceId face_id = node_to_face_matrix[n_id][i_face];
+ for (size_t i_node = 0; i_node < face_to_node_matrix[face_id].size(); ++i_node) {
+ NodeId node_id = face_to_node_matrix[face_id][i_node];
+ if (node_id == n_id) {
+ continue;
+ } else if ((is_fixed[node_id])) {
+ Afree(local_free_id, local_free_id) += 1;
+ Afixed(local_free_id, id_fixed[node_id]) -= 1;
+ } else {
+ Afree(local_free_id, local_free_id) += 1;
+ Afree(local_free_id, id_free[node_id]) -= 1;
+ }
}
}
}
@@ -239,6 +265,49 @@ class ImplicitMeshSmootherHandler::ImplicitMeshSmoother
}
}
+ void
+ _correct_sym(const NodeValue<const Rd>& old_xr, NodeValue<Rd>& new_xr) const
+ {
+ for (auto&& boundary_condition : m_boundary_condition_list) {
+ std::visit(
+ [&](auto&& bc) {
+ using BCType = std::decay_t<decltype(bc)>;
+ if constexpr (std::is_same_v<BCType, SymmetryBoundaryCondition>) {
+ const Rd& n = bc.outgoingNormal();
+
+ const Rdxd I = identity;
+ const Rdxd nxn = tensorProduct(n, n);
+ const Rdxd P = I - nxn;
+
+ const Array<const NodeId>& node_list = bc.nodeList();
+ parallel_for(
+ node_list.size(), PUGS_LAMBDA(const size_t i_node) {
+ const NodeId node_id = node_list[i_node];
+ Rd depl = new_xr[node_id] - old_xr[node_id];
+ new_xr[node_id] = old_xr[node_id] + P * depl;
+ });
+ // } else if constexpr (std::is_same_v<BCType, AxisBoundaryCondition>) {
+ // if constexpr (Dimension > 1) {
+ // const Rd& t = bc.direction();
+
+ // const Rdxd txt = tensorProduct(t, t);
+
+ // const Array<const NodeId>& node_list = bc.nodeList();
+ // parallel_for(
+ // node_list.size(), PUGS_LAMBDA(const size_t i_node) {
+ // const NodeId node_id = node_list[i_node];
+
+ // shift[node_id] = txt * shift[node_id];
+ // });
+ // } else {
+ // throw UnexpectedError("AxisBoundaryCondition make no sense in dimension 1");
+ // }
+ }
+ },
+ boundary_condition);
+ }
+ }
+
NodeValue<Rd>
_getPosition(NodeValue<const bool> is_displaced) const
{
@@ -247,7 +316,9 @@ class ImplicitMeshSmootherHandler::ImplicitMeshSmoother
NodeValue<Rd> pos_r{connectivity};
NodeValue<bool> is_fixed{connectivity};
+ NodeValue<int> is_symmetric{connectivity};
is_fixed.fill(false);
+ is_symmetric.fill(0);
parallel_for(
m_given_mesh.numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
pos_r[node_id] = given_xr[node_id];
@@ -255,24 +326,26 @@ class ImplicitMeshSmootherHandler::ImplicitMeshSmoother
is_fixed[node_id] = true;
}
});
- _browseBC(is_fixed);
+ _browseBC(is_fixed, is_symmetric);
+ // std::cout << "is_fixed " << is_fixed << "\n";
+ // std::cout << "is_symmetric " << is_symmetric << "\n";
NodeValue<int> node_dof_id{connectivity};
size_t nb_free, nb_fixed;
NodeValue<int> id_free{connectivity};
NodeValue<int> id_fixed{connectivity};
_computeMatrixSize(is_fixed, id_fixed, id_free, nb_free, nb_fixed);
+ std::cout << " nb_free " << nb_free << " nb_fixed " << nb_fixed << "\n";
Array<int> non_zeros_free{nb_free};
Array<int> non_zeros_fixed{nb_free};
Array<NodeId> gid_node_free{nb_free};
Array<NodeId> gid_node_fixed{nb_fixed};
_findLocalIds(is_fixed, gid_node_free, gid_node_fixed, non_zeros_free, non_zeros_fixed);
- std::cout << " nb_free " << nb_free << " nb_fixed " << nb_fixed << "\n";
CRSMatrixDescriptor<double> Afree(nb_free, nb_free, non_zeros_free);
CRSMatrixDescriptor<double> Afixed(nb_free, nb_fixed, non_zeros_fixed);
LinearSolver solver;
- _fillMatrix(Afree, Afixed, is_fixed, id_free, id_fixed);
+ _fillMatrix(Afree, Afixed, is_fixed, is_symmetric, id_free, id_fixed);
Vector<double> F{nb_fixed};
CRSMatrix Mfree{Afree.getCRSMatrix()};
CRSMatrix Mfixed{Afixed.getCRSMatrix()};
@@ -283,6 +356,10 @@ class ImplicitMeshSmootherHandler::ImplicitMeshSmoother
Vector<double> X{nb_free};
Vector<double> b{nb_free};
b = Mfixed * F;
+ // std::cout << " Mfixed " << Mfixed << "\n";
+ // std::cout << " F " << F << "\n";
+ // std::cout << " b " << b << "\n";
+ // std::cout << " Mfree " << Mfree << "\n";
solver.solveLocalSystem(Mfree, X, b);
parallel_for(
m_given_mesh.numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
@@ -291,6 +368,7 @@ class ImplicitMeshSmootherHandler::ImplicitMeshSmoother
}
});
}
+ _correct_sym(m_given_mesh.xr(), pos_r);
// synchronize(pos_r);
return pos_r;
}