@@ -43,16 +43,15 @@ class HyperElasticProblem
4343{
4444public:
4545 // / Constructor
46- HyperElasticProblem (
47- std::shared_ptr<fem::Form<T>> L, std::shared_ptr<fem::Form<T>> J,
48- std::vector<std::shared_ptr<const fem::DirichletBC<T>>> bcs)
49- : _l(L), _j(J), _bcs(bcs),
50- _b (L->function_spaces ()[0]->dofmap()->index_map,
51- L->function_spaces()[0]->dofmap()->index_map_bs()),
52- _matA(la::petsc::Matrix(fem::petsc::create_matrix(*J, " aij" ), false))
46+ HyperElasticProblem (fem::Form<T>& L, fem::Form<T>& J,
47+ const std::vector<fem::DirichletBC<T>>& bcs)
48+ : _l(L), _j(J), _bcs(bcs.begin(), bcs.end()),
49+ _b (L.function_spaces()[0]->dofmap()->index_map,
50+ L.function_spaces()[0]->dofmap()->index_map_bs()),
51+ _matA(la::petsc::Matrix(fem::petsc::create_matrix(J, " aij" ), false))
5352 {
54- auto map = L-> function_spaces ()[0 ]->dofmap ()->index_map ;
55- const int bs = L-> function_spaces ()[0 ]->dofmap ()->index_map_bs ();
53+ auto map = L. function_spaces ()[0 ]->dofmap ()->index_map ;
54+ const int bs = L. function_spaces ()[0 ]->dofmap ()->index_map_bs ();
5655 std::int32_t size_local = bs * map->size_local ();
5756
5857 std::vector<PetscInt> ghosts (map->ghosts ().begin (), map->ghosts ().end ());
@@ -88,7 +87,7 @@ class HyperElasticProblem
8887 // Assemble b and update ghosts
8988 std::span b (_b.mutable_array ());
9089 std::ranges::fill (b, 0 );
91- fem::assemble_vector<T>(b, * _l);
90+ fem::assemble_vector<T>(b, _l);
9291 VecGhostUpdateBegin (_b_petsc, ADD_VALUES , SCATTER_REVERSE );
9392 VecGhostUpdateEnd (_b_petsc, ADD_VALUES , SCATTER_REVERSE );
9493
@@ -100,7 +99,7 @@ class HyperElasticProblem
10099 const T* _x = nullptr ;
101100 VecGetArrayRead (x_local, &_x);
102101 std::ranges::for_each (_bcs, [b, x = std::span (_x, n)](auto & bc)
103- { bc-> set (b, x, -1 ); });
102+ { bc. get (). set (b, x, -1 ); });
104103 VecRestoreArrayRead (x_local, &_x);
105104 };
106105 }
@@ -111,12 +110,12 @@ class HyperElasticProblem
111110 return [&](const Vec, Mat A)
112111 {
113112 MatZeroEntries (A);
114- fem::assemble_matrix (la::petsc::Matrix::set_block_fn (A, ADD_VALUES ), * _j,
113+ fem::assemble_matrix (la::petsc::Matrix::set_block_fn (A, ADD_VALUES ), _j,
115114 _bcs);
116115 MatAssemblyBegin (A, MAT_FLUSH_ASSEMBLY );
117116 MatAssemblyEnd (A, MAT_FLUSH_ASSEMBLY );
118117 fem::set_diagonal (la::petsc::Matrix::set_fn (A, INSERT_VALUES ),
119- *_j-> function_spaces ()[0 ], _bcs);
118+ *_j. function_spaces ()[0 ], _bcs);
120119 MatAssemblyBegin (A, MAT_FINAL_ASSEMBLY );
121120 MatAssemblyEnd (A, MAT_FINAL_ASSEMBLY );
122121 };
@@ -129,8 +128,9 @@ class HyperElasticProblem
129128 Mat matrix () { return _matA.mat (); }
130129
131130private:
132- std::shared_ptr<fem::Form<T>> _l, _j;
133- std::vector<std::shared_ptr<const fem::DirichletBC<T>>> _bcs;
131+ fem::Form<T>& _l;
132+ fem::Form<T>& _j;
133+ std::vector<std::reference_wrapper<const fem::DirichletBC<T>>> _bcs;
134134 la::Vector<T> _b;
135135 Vec _b_petsc = nullptr ;
136136 la::petsc::Matrix _matA;
@@ -166,20 +166,22 @@ int main(int argc, char* argv[])
166166 basix::element::lagrange_variant::unset,
167167 basix::element::dpc_variant::unset, false );
168168
169- auto V = std::make_shared<fem::FunctionSpace<U>>(
170- fem::create_functionspace (mesh, element, {3 }));
169+ auto V
170+ = std::make_shared<fem::FunctionSpace<U>>(fem::create_functionspace<U>(
171+ mesh, std::make_shared<fem::FiniteElement<U>>(
172+ element, std::vector<std::size_t >{3 })));
171173
172174 auto B = std::make_shared<fem::Constant<T>>(std::vector<T>{0 , 0 , 0 });
173175 auto traction = std::make_shared<fem::Constant<T>>(std::vector<T>{0 , 0 , 0 });
174176
175177 // Define solution function
176178 auto u = std::make_shared<fem::Function<T>>(V);
177- auto a = std::make_shared< fem::Form<T>>(
178- fem::create_form<T>(*form_hyperelasticity_J_form, {V, V}, {{" u" , u}},
179- {{" B" , B}, {" T" , traction}}, {}, {}) );
180- auto L = std::make_shared< fem::Form<T>>(
181- fem::create_form<T>(*form_hyperelasticity_F_form, {V}, {{" u" , u}},
182- {{" B" , B}, {" T" , traction}}, {}, {}) );
179+ fem::Form<T> a
180+ = fem::create_form<T>(*form_hyperelasticity_J_form, {V, V}, {{" u" , u}},
181+ {{" B" , B}, {" T" , traction}}, {}, {});
182+ fem::Form<T> L
183+ = fem::create_form<T>(*form_hyperelasticity_F_form, {V}, {{" u" , u}},
184+ {{" B" , B}, {" T" , traction}}, {}, {});
183185
184186 auto u_rotation = std::make_shared<fem::Function<T>>(V);
185187 u_rotation->interpolate (
@@ -243,10 +245,9 @@ int main(int argc, char* argv[])
243245 }
244246 return marker;
245247 });
246- std::vector bcs = {
247- std::make_shared<const fem::DirichletBC<T>>(std::vector<T>{0 , 0 , 0 },
248- bdofs_left, V),
249- std::make_shared<const fem::DirichletBC<T>>(u_rotation, bdofs_right)};
248+ std::vector bcs
249+ = {fem::DirichletBC<T>(std::vector<T>{0 , 0 , 0 }, bdofs_left, V),
250+ fem::DirichletBC<T>(u_rotation, bdofs_right)};
250251
251252 HyperElasticProblem problem (L, a, bcs);
252253 nls::petsc::NewtonSolver newton_solver (mesh->comm ());
@@ -262,6 +263,9 @@ int main(int argc, char* argv[])
262263
263264 // Compute Cauchy stress. Construct appropriate Basix element for
264265 // stress.
266+ fem::Expression sigma_expression = fem::create_expression<T, U>(
267+ *expression_hyperelasticity_sigma, {{" u" , u}}, {});
268+
265269 constexpr auto family = basix::element::family::P;
266270 auto cell_type
267271 = mesh::cell_type_to_basix_type (mesh->topology ()->cell_type ());
@@ -270,12 +274,12 @@ int main(int argc, char* argv[])
270274 basix::FiniteElement S_element = basix::create_element<U>(
271275 family, cell_type, k, basix::element::lagrange_variant::unset,
272276 basix::element::dpc_variant::unset, discontinuous);
273- auto S = std::make_shared<fem::FunctionSpace<U>>( fem::create_functionspace (
274- mesh, S_element, std::vector<std:: size_t >{ 3 , 3 }));
275- auto sigma_expression = fem::create_expression<T, U >(
276- *expression_hyperelasticity_sigma, {{ " u " , u}}, {} );
277+ auto S
278+ = std::make_shared<fem::FunctionSpace<U>>(fem::create_functionspace<U>(
279+ mesh, std::make_shared< fem::FiniteElement<U> >(
280+ S_element, std::vector<std:: size_t >{ 3 , 3 })) );
277281
278- auto sigma = fem::Function<T>(S);
282+ fem::Function<T> sigma (S);
279283 sigma.name = " cauchy_stress" ;
280284 sigma.interpolate (sigma_expression);
281285
0 commit comments