50 if (formulation ==
"Helmholtz")
51 return std::make_shared<Helmholtz>();
52 else if (formulation ==
"Laplacian")
53 return std::make_shared<Laplacian>();
54 else if (formulation ==
"Bilaplacian")
55 return std::make_shared<BilaplacianMain>();
56 else if (formulation ==
"BilaplacianAux")
57 return std::make_shared<BilaplacianAux>();
59 else if (formulation ==
"LinearElasticity")
60 return std::make_shared<LinearElasticity>();
61 else if (formulation ==
"HookeLinearElasticity")
62 return std::make_shared<HookeLinearElasticity>();
63 else if (formulation ==
"IncompressibleLinearElasticity")
64 return std::make_shared<IncompressibleLinearElasticityDispacement>();
65 else if (formulation ==
"IncompressibleLinearElasticityPressure")
66 return std::make_shared<IncompressibleLinearElasticityPressure>();
68 else if (formulation ==
"SaintVenant")
69 return std::make_shared<SaintVenantElasticity>();
70 else if (formulation ==
"NeoHookean")
71 return std::make_shared<NeoHookeanElasticity>();
72 else if (formulation ==
"MooneyRivlin")
73 return std::make_shared<MooneyRivlinElasticity>();
74 else if (formulation ==
"MooneyRivlin3Param")
75 return std::make_shared<MooneyRivlin3ParamElasticity>();
76 else if (formulation ==
"MooneyRivlin3ParamSymbolic")
77 return std::make_shared<MooneyRivlin3ParamSymbolic>();
78 else if (formulation ==
"MultiModels")
79 return std::make_shared<MultiModel>();
80 else if (formulation ==
"UnconstrainedOgden")
81 return std::make_shared<UnconstrainedOgdenElasticity>();
82 else if (formulation ==
"IncompressibleOgden")
83 return std::make_shared<IncompressibleOgdenElasticity>();
85 else if (formulation ==
"Stokes")
86 return std::make_shared<StokesVelocity>();
87 else if (formulation ==
"StokesPressure")
88 return std::make_shared<StokesPressure>();
89 else if (formulation ==
"NavierStokes")
90 return std::make_shared<NavierStokesVelocity>();
91 else if (formulation ==
"OperatorSplitting")
92 return std::make_shared<OperatorSplitting>();
94 else if (formulation ==
"AMIPS")
95 return std::make_shared<AMIPSEnergy>();
113 const int n_bases,
const int n_pressure_bases,
const int problem_dim,
const bool add_average,
117 assert(velocity_stiffness.rows() == velocity_stiffness.cols());
118 assert(velocity_stiffness.rows() == n_bases * problem_dim);
120 assert(mixed_stiffness.size() == 0 || mixed_stiffness.rows() == n_bases * problem_dim);
121 assert(mixed_stiffness.size() == 0 || mixed_stiffness.cols() == n_pressure_bases);
123 assert(pressure_stiffness.size() == 0 || pressure_stiffness.rows() == n_pressure_bases);
124 assert(pressure_stiffness.size() == 0 || pressure_stiffness.cols() == n_pressure_bases);
126 const int avg_offset = add_average ? 1 : 0;
128 std::vector<Eigen::Triplet<double>> blocks;
129 blocks.reserve(velocity_stiffness.nonZeros() + 2 * mixed_stiffness.nonZeros() + pressure_stiffness.nonZeros() + 2 * avg_offset * velocity_stiffness.rows());
131 for (
int k = 0; k < velocity_stiffness.outerSize(); ++k)
133 for (StiffnessMatrix::InnerIterator it(velocity_stiffness, k); it; ++it)
135 blocks.emplace_back(it.row(), it.col(), it.value());
139 for (
int k = 0; k < mixed_stiffness.outerSize(); ++k)
141 for (StiffnessMatrix::InnerIterator it(mixed_stiffness, k); it; ++it)
143 blocks.emplace_back(it.row(), n_bases * problem_dim + it.col(), it.value());
144 blocks.emplace_back(it.col() + n_bases * problem_dim, it.row(), it.value());
148 for (
int k = 0; k < pressure_stiffness.outerSize(); ++k)
150 for (StiffnessMatrix::InnerIterator it(pressure_stiffness, k); it; ++it)
152 blocks.emplace_back(n_bases * problem_dim + it.row(), n_bases * problem_dim + it.col(), it.value());
158 const double val = 1.0 / n_pressure_bases;
159 for (
int i = 0; i < n_pressure_bases; ++i)
161 blocks.emplace_back(n_bases * problem_dim + i, n_bases * problem_dim + n_pressure_bases,
val);
162 blocks.emplace_back(n_bases * problem_dim + n_pressure_bases, n_bases * problem_dim + i,
val);
166 stiffness.resize(n_bases * problem_dim + n_pressure_bases + avg_offset, n_bases * problem_dim + n_pressure_bases + avg_offset);
167 stiffness.setFromTriplets(blocks.begin(), blocks.end());
168 stiffness.makeCompressed();
180 if (assembler ==
"Mass")
184 return std::max(basis_degree * 2, 1);
186 return basis_degree * 2 + 1;
188 else if (assembler ==
"NavierStokes")
191 return std::max((basis_degree - 1) + basis_degree, 1);
193 return std::max(basis_degree * 2, 1);
195 return basis_degree * 2 + 1;
202 return std::max((basis_degree - 1) * 2, 1);
212 return std::max(basis_degree * 2, 1);
215 return (basis_degree - 1) * 2 + 1;
static int quadrature_order(const std::string &assembler, const int basis_degree, const BasisType &b_type, const int dim)
utility for retrieving the needed quadrature order to precisely integrate the given form on the given...
static void merge_mixed_matrices(const int n_bases, const int n_pressure_bases, const int problem_dim, const bool add_average, const StiffnessMatrix &velocity_stiffness, const StiffnessMatrix &mixed_stiffness, const StiffnessMatrix &pressure_stiffness, StiffnessMatrix &stiffness)
utility to merge 3 blocks of mixed matrices, A=velocity_stiffness, B=mixed_stiffness,...