9 Eigen::MatrixXd extract_vertices(
const std::shared_ptr<polyfem::mesh::MeshNodes> &mesh_nodes)
12 V.setZero(mesh_nodes->n_nodes(), mesh_nodes->node_position(0).size());
14 for (
int i = 0; i <
V.rows(); i++)
15 V.row(i) = mesh_nodes->node_position(i);
26 const std::vector<basis::ElementBases> &bases,
27 const std::shared_ptr<mesh::MeshNodes> &mesh_nodes,
28 const Eigen::MatrixXd &affine_matrix,
29 const double tol): affine_matrix_(affine_matrix)
31 const int dim = affine_matrix.rows();
34 Eigen::MatrixXd
V = extract_vertices(mesh_nodes) * affine_matrix.inverse().transpose();
35 auto boundary_nodes = mesh_nodes->boundary_nodes();
39 const double bbox_size = (max - min).maxCoeff();
40 for (
int d = 0; d <
dim; d++)
42 const double eps = tol * bbox_size;
52 for (
int d = 0; d <
dim; d++)
54 std::set<int> dependent_ids, target_ids;
55 for (
int bid : boundary_nodes)
57 if (
V(bid, d) > max(d) - eps)
58 target_ids.insert(bid);
59 else if (
V(bid, d) < min(d) + eps)
60 dependent_ids.insert(bid);
63 for (
int i : dependent_ids)
65 bool found_target =
false;
66 for (
int j : target_ids)
69 projected_diff(d) = 0;
70 if (projected_diff.norm() < eps)
96 for (
int d = 0; d <
dim; d++)
108 int independent_dof = 0;
121 for (
int i = 0; i < old_size; i++)
131 auto extended_index_map = [&](
int id){
138 StiffnessMatrix A_periodic(extended_index_map(A.rows()), extended_index_map(A.cols()));
139 std::vector<Eigen::Triplet<double>>
entries;
141 for (
int k = 0; k < A.outerSize(); k++)
143 for (StiffnessMatrix::InnerIterator it(A,k); it; ++it)
145 entries.emplace_back(extended_index_map(it.row()), extended_index_map(it.col()), it.value());
150 std::swap(A_periodic, A);
152 return independent_dof;
159 auto extended_index_map = [&](
int id){
167 Eigen::MatrixXd b_periodic;
168 b_periodic.setZero(extended_index_map(b.rows()), b.cols());
170 for (
int k = 0; k < b.rows(); k++)
171 b_periodic.row(extended_index_map(k)) += b.row(k);
173 for (
int k = 0; k < b.rows(); k++)
174 b_periodic.row(extended_index_map(k)) = b.row(k);
180 std::vector<int> boundary_nodes_reduced = boundary_nodes;
183 for (
int i = 0; i < boundary_nodes_reduced.size(); i++)
186 std::sort(boundary_nodes_reduced.begin(), boundary_nodes_reduced.end());
187 auto it = std::unique(boundary_nodes_reduced.begin(), boundary_nodes_reduced.end());
188 boundary_nodes_reduced.resize(std::distance(boundary_nodes_reduced.begin(), it));
190 return boundary_nodes_reduced;
197 auto extended_index_map = [&](
int id){
204 Eigen::MatrixXd x_full;
205 x_full.resize(ndofs, x_periodic.cols());
206 for (
int i = 0; i < x_full.rows(); i++)
207 x_full.row(i) = x_periodic.row(extended_index_map(i));
std::vector< Eigen::Triplet< double > > entries
Eigen::MatrixXd affine_matrix_
Eigen::VectorXi dependent_map_
std::vector< std::array< int, 2 > > constraint_list_
Eigen::VectorXi periodic_mask_
Eigen::VectorXi full_to_periodic_map_
bool has_periodic_bc() const
bool all_direction_periodic() const
Eigen::MatrixXd periodic_to_full(const int ndofs, const Eigen::MatrixXd &x_periodic) const
PeriodicBoundary(const bool is_scalar, const int n_bases, const std::vector< basis::ElementBases > &bases, const std::shared_ptr< mesh::MeshNodes > &mesh_nodes, const Eigen::MatrixXd &affine_matrix, const double tol)
Eigen::VectorXi index_map_
int full_to_periodic(StiffnessMatrix &A) const
Eigen::Matrix< double, 1, Eigen::Dynamic, Eigen::RowMajor, 1, 3 > RowVectorNd
void log_and_throw_error(const std::string &msg)
Eigen::SparseMatrix< double, Eigen::ColMajor > StiffnessMatrix