14 const std::vector<std::shared_ptr<Form>> &forms,
15 const std::vector<std::shared_ptr<AugmentedLagrangianForm>> &penalty_forms,
16 const bool solve_symmetric_macro_strain,
17 const std::shared_ptr<polysolve::linear::Solver> &solver)
18 :
NLProblem(full_size, state.periodic_bc, t, forms, penalty_forms, solver),
20 only_symmetric(solve_symmetric_macro_strain),
21 macro_strain_constraint_(macro_strain_constraint)
33 for (
int i = 0, idx = 0; i < dim; i++)
35 for (
int j = i; j < dim; j++)
60 Eigen::VectorXd reduced(dof1 + dof2);
72 assert(reduced.size() == dof1 + dof2);
76 Eigen::VectorXd extended(full.size() + disp_grad.size());
77 extended << full, disp_grad;
88 Eigen::VectorXd grad(dof1 + dof2);
102 term.segment(i * dim, dim) += disp_grad.transpose() * adjoint_full.segment(i * dim, dim);
140 Eigen::VectorXd grad_extended;
151 Eigen::MatrixXd A12, A22;
153 Eigen::MatrixXd tmp = Eigen::MatrixXd(extended.rightCols(dim * dim));
158 std::vector<Eigen::Triplet<double>>
entries;
159 entries.reserve(extended.nonZeros() + A12.size() * 2 + A22.size());
162 for (StiffnessMatrix::InnerIterator it(extended, k); it; ++it)
164 entries.emplace_back(it.row(), it.col(), it.value());
166 for (
int i = 0; i < A12.rows(); i++)
168 for (
int j = 0; j < A12.cols(); j++)
175 for (
int i = 0; i < A22.rows(); i++)
176 for (
int j = 0; j < A22.cols(); j++)
183 reduced.resize(0, 0);
188 assert(dof1 ==
Q2_.cols());
193 for (
int k = 0; k <
Q2_.cols(); ++k)
194 for (StiffnessMatrix::InnerIterator it(
Q2_, k); it; ++it)
196 entries.emplace_back(it.row(), it.col(), it.value());
199 for (
int k = 0; k < dof2; k++)
201 entries.emplace_back(
Q2_.rows() + k, dof1 + k, 1);
207 reduced = Q2_extended.transpose() * reduced * Q2_extended;
209 reduced.prune([](
const Eigen::Index &row,
const Eigen::Index &col,
const Scalar &
value) {
210 return std::abs(
value) > 1e-10;
231 THessian hess_extended;
244 Eigen::VectorXd fixed_mask;
245 fixed_mask.setZero(dim * dim);
246 fixed_mask(fixed_entry.array()).setOnes();
250 for (
int i = 0; i < fixed_mask.size(); i++)
251 if (abs(fixed_mask(i)) > 1e-8)
256 for (
int i = 0, j = 0; i <
fixed_mask_.size(); i++)
269 Eigen::MatrixXd A12 =
hessian * tmp.transpose();
270 Eigen::MatrixXd A22 = tmp * A12;
272 std::vector<Eigen::Triplet<double>>
entries;
273 entries.reserve(
hessian.nonZeros() + A12.size() * 2 + A22.size());
275 for (
int k = 0; k <
hessian.outerSize(); ++k)
276 for (StiffnessMatrix::InnerIterator it(
hessian, k); it; ++it)
277 entries.emplace_back(it.row(), it.col(), it.value());
279 for (
int i = 0; i < A12.rows(); i++)
280 for (
int j = 0; j < A12.cols(); j++)
286 for (
int i = 0; i < A22.rows(); i++)
287 for (
int j = 0; j < A22.cols(); j++)
296 assert(dof1 ==
Q2_.cols());
301 for (
int k = 0; k <
Q2_.cols(); ++k)
302 for (StiffnessMatrix::InnerIterator it(
Q2_, k); it; ++it)
304 entries.emplace_back(it.row(), it.col(), it.value());
307 for (
int k = 0; k < dof2; k++)
309 entries.emplace_back(
Q2_.rows() + k, dof1 + k, 1);
317 hessian.prune([](
const Eigen::Index &row,
const Eigen::Index &col,
const Scalar &
value) {
318 return std::abs(
value) > 1e-10;
336 reduced.setZero(dof1 + dof2);
350 reduced.setZero(dof1 + dof2);
384 mid(i) = fixed_values(i);
433 form->post_step(polysolve::nonlinear::PostStepData(
472 for (
int i = 0; i < X.rows(); i++)
473 for (
int j = 0; j < dim; j++)
474 for (
int k = 0; k < dim; k++)
475 jac(j * dim + k, i * dim + j) = X(i, k);
std::vector< Eigen::Triplet< double > > entries
main class that contains the polyfem solver and all its state
int n_bases
number of bases
std::unique_ptr< mesh::Mesh > mesh
current mesh, it can be a Mesh2D or Mesh3D
std::shared_ptr< polyfem::mesh::MeshNodes > mesh_nodes
Mapping from input nodes to FE nodes.
StiffnessMatrix basis_nodes_to_gbasis_nodes
Eigen::MatrixXd eval(const double t) const
static Eigen::MatrixXd generate_linear_field(const int n_bases, const std::shared_ptr< mesh::MeshNodes > mesh_nodes, const Eigen::MatrixXd &grad)
static Eigen::MatrixXd get_bases_position(const int n_bases, const std::shared_ptr< mesh::MeshNodes > mesh_nodes)
virtual void hessian(const TVector &x, THessian &hessian) override
virtual double value(const TVector &x) override
virtual void init(const TVector &x0) override
virtual void gradient(const TVector &x, TVector &gradv) override
Eigen::MatrixXd reduced_to_disp_grad(const TVector &reduced, bool homogeneous=false) const
Eigen::VectorXi fixed_mask_
void gradient(const TVector &x, TVector &gradv) override
void init(const TVector &x0) override
void post_step(const polysolve::nonlinear::PostStepData &data) override
Eigen::MatrixXd macro_mid_to_reduced_
int macro_reduced_size() const
void set_fixed_entry(const Eigen::VectorXi &fixed_entry)
TVector macro_full_to_reduced(const TVector &full) const
bool is_step_valid(const TVector &x0, const TVector &x1) override
Eigen::MatrixXd macro_mid_to_full_
const assembler::MacroStrainValue & macro_strain_constraint_
void line_search_begin(const TVector &x0, const TVector &x1) override
void update_quantities(const double t, const TVector &x) override
Eigen::MatrixXd macro_full_to_mid_
TVector extended_to_reduced_grad(const TVector &extended) const
NLHomoProblem(const int full_size, const assembler::MacroStrainValue ¯o_strain_constraint, const State &state, const double t, const std::vector< std::shared_ptr< Form > > &forms, const std::vector< std::shared_ptr< AugmentedLagrangianForm > > &penalty_forms, const bool solve_symmetric_macro_strain, const std::shared_ptr< polysolve::linear::Solver > &solver)
void hessian(const TVector &x, THessian &hessian) override
void init_lagging(const TVector &x) override
TVector reduced_to_extended(const TVector &reduced, bool homogeneous=false) const
TVector reduced_to_full(const TVector &reduced) const
TVector extended_to_reduced(const TVector &extended) const
std::vector< std::shared_ptr< Form > > homo_forms
Eigen::MatrixXd constraint_grad() const
TVector reduced_to_full_shape_derivative(const Eigen::MatrixXd &disp_grad, const TVector &adjoint_full) const
Eigen::MatrixXd macro_full_to_reduced_grad(const Eigen::MatrixXd &full) const
double value(const TVector &x) override
TVector full_to_reduced(const TVector &full, const Eigen::MatrixXd &disp_grad) const
TVector full_to_reduced_grad(const TVector &full) const override
void full_hessian_to_reduced_hessian(THessian &hessian) const
const bool only_symmetric
void solution_changed(const TVector &new_x) override
double max_step_size(const TVector &x0, const TVector &x1) override
void extended_hessian_to_reduced_hessian(const THessian &extended, THessian &reduced) const
bool is_step_collision_free(const TVector &x0, const TVector &x1) override
void update_lagging(const TVector &x, const int iter_num) override
TVector macro_reduced_to_full(const TVector &reduced, bool homogeneous=false) const
const int full_size_
Size of the full problem.
void line_search_begin(const TVector &x0, const TVector &x1) override
StiffnessMatrix Q2_
Q2 block of the QR decomposition of the constraints matrix.
virtual bool is_step_valid(const TVector &x0, const TVector &x1) override
virtual void post_step(const polysolve::nonlinear::PostStepData &data) override
virtual TVector full_to_reduced_grad(const TVector &full) const
TVector full_to_reduced(const TVector &full) const
virtual void update_quantities(const double t, const TVector &x)
void init_lagging(const TVector &x) override
virtual bool is_step_collision_free(const TVector &x0, const TVector &x1) override
TVector reduced_to_full(const TVector &reduced) const
void update_lagging(const TVector &x, const int iter_num) override
virtual double max_step_size(const TVector &x0, const TVector &x1) override
void solution_changed(const TVector &new_x) override
std::shared_ptr< FullNLProblem > penalty_problem_
Eigen::MatrixXd unflatten(const Eigen::VectorXd &x, int dim)
Unflatten rowwises, so every dim elements in x become a row.
Eigen::VectorXd flatten(const Eigen::MatrixXd &X)
Flatten rowwises.
void log_and_throw_error(const std::string &msg)
Eigen::SparseMatrix< double, Eigen::ColMajor > StiffnessMatrix