Loading [MathJax]/extensions/tex2jax.js
PolyFEM
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
HookeLinearElasticity.cpp
Go to the documentation of this file.
2
4
5namespace polyfem::assembler
6{
7 using namespace basis;
8 namespace
9 {
10 template <class Matrix>
11 Matrix strain_from_disp_grad(const Matrix &disp_grad)
12 {
13 Matrix mat = (disp_grad + disp_grad.transpose());
14
15 for (int i = 0; i < mat.size(); ++i)
16 mat(i) *= 0.5;
17
18 return mat;
19 }
20
21 template <int dim>
22 Eigen::Matrix<double, dim, dim> strain(const Eigen::MatrixXd &grad, const Eigen::MatrixXd &jac_it, int k, int coo)
23 {
24 Eigen::Matrix<double, dim, dim> jac;
25 jac.setZero();
26 jac.row(coo) = grad.row(k);
27 jac = jac * jac_it;
28
29 return strain_from_disp_grad(jac);
30 }
31
32 template <typename T, unsigned long N>
33 T stress(const ElasticityTensor &elasticity_tensor, const std::array<T, N> &strain, const int j)
34 {
35 T res = elasticity_tensor(j, 0) * strain[0];
36
37 for (unsigned long k = 1; k < N; ++k)
38 res += elasticity_tensor(j, k) * strain[k];
39
40 return res;
41 }
42 } // namespace
43
47
48 void HookeLinearElasticity::add_multimaterial(const int index, const json &params, const Units &units)
49 {
50 assert(size() == 2 || size() == 3);
51
52 if (!params.contains("elasticity_tensor") || params["elasticity_tensor"].empty())
53 {
54 if (params.count("young"))
55 {
56 if (params["young"].is_number() && params["nu"].is_number())
57 elasticity_tensor_.set_from_young_poisson(params["young"], params["nu"], units.stress());
58 }
59 else if (params.count("E"))
60 {
61 if (params["E"].is_number() && params["nu"].is_number())
62 elasticity_tensor_.set_from_young_poisson(params["E"], params["nu"], units.stress());
63 }
64 else if (params.count("lambda"))
65 {
66 if (params["lambda"].is_number() && params["mu"].is_number())
67 elasticity_tensor_.set_from_lambda_mu(params["lambda"], params["mu"], units.stress());
68 }
69 }
70 else
71 {
72 std::vector<double> entries = params["elasticity_tensor"];
74 }
75 }
76
82
83 Eigen::Matrix<double, Eigen::Dynamic, 1, 0, 9, 1>
85 {
86 const Eigen::MatrixXd &gradi = data.vals.basis_values[data.i].grad;
87 const Eigen::MatrixXd &gradj = data.vals.basis_values[data.j].grad;
88
89 // (C : gradi) : gradj
90 Eigen::Matrix<double, Eigen::Dynamic, 1, 0, 9, 1> res(size() * size());
91 res.setZero();
92 assert(gradi.cols() == size());
93 assert(gradj.cols() == size());
94 assert(size_t(gradi.rows()) == data.vals.jac_it.size());
95
96 for (long k = 0; k < gradi.rows(); ++k)
97 {
98 Eigen::Matrix<double, Eigen::Dynamic, 1, 0, 9, 1> res_k(size() * size());
99
100 if (size() == 2)
101 {
102 const Eigen::Matrix2d eps_x_i = strain<2>(gradi, data.vals.jac_it[k], k, 0);
103 const Eigen::Matrix2d eps_y_i = strain<2>(gradi, data.vals.jac_it[k], k, 1);
104
105 const Eigen::Matrix2d eps_x_j = strain<2>(gradj, data.vals.jac_it[k], k, 0);
106 const Eigen::Matrix2d eps_y_j = strain<2>(gradj, data.vals.jac_it[k], k, 1);
107
108 std::array<double, 3> e_x, e_y;
109 e_x[0] = eps_x_i(0, 0);
110 e_x[1] = eps_x_i(1, 1);
111 e_x[2] = 2 * eps_x_i(0, 1);
112
113 e_y[0] = eps_y_i(0, 0);
114 e_y[1] = eps_y_i(1, 1);
115 e_y[2] = 2 * eps_y_i(0, 1);
116
117 Eigen::Matrix2d sigma_x;
118 sigma_x << elasticity_tensor_.compute_stress<3>(e_x, 0), elasticity_tensor_.compute_stress<3>(e_x, 2),
120
121 Eigen::Matrix2d sigma_y;
122 sigma_y << elasticity_tensor_.compute_stress<3>(e_y, 0), elasticity_tensor_.compute_stress<3>(e_y, 2),
124
125 res_k(0) = (sigma_x * eps_x_j).trace();
126 res_k(2) = (sigma_x * eps_y_j).trace();
127
128 res_k(1) = (sigma_y * eps_x_j).trace();
129 res_k(3) = (sigma_y * eps_y_j).trace();
130 }
131 else
132 {
133 const Eigen::Matrix3d eps_x_i = strain<3>(gradi, data.vals.jac_it[k], k, 0);
134 const Eigen::Matrix3d eps_y_i = strain<3>(gradi, data.vals.jac_it[k], k, 1);
135 const Eigen::Matrix3d eps_z_i = strain<3>(gradi, data.vals.jac_it[k], k, 2);
136
137 const Eigen::Matrix3d eps_x_j = strain<3>(gradj, data.vals.jac_it[k], k, 0);
138 const Eigen::Matrix3d eps_y_j = strain<3>(gradj, data.vals.jac_it[k], k, 1);
139 const Eigen::Matrix3d eps_z_j = strain<3>(gradj, data.vals.jac_it[k], k, 2);
140
141 std::array<double, 6> e_x, e_y, e_z;
142 e_x[0] = eps_x_i(0, 0);
143 e_x[1] = eps_x_i(1, 1);
144 e_x[2] = eps_x_i(2, 2);
145 e_x[3] = 2 * eps_x_i(1, 2);
146 e_x[4] = 2 * eps_x_i(0, 2);
147 e_x[5] = 2 * eps_x_i(0, 1);
148
149 e_y[0] = eps_y_i(0, 0);
150 e_y[1] = eps_y_i(1, 1);
151 e_y[2] = eps_y_i(2, 2);
152 e_y[3] = 2 * eps_y_i(1, 2);
153 e_y[4] = 2 * eps_y_i(0, 2);
154 e_y[5] = 2 * eps_y_i(0, 1);
155
156 e_z[0] = eps_z_i(0, 0);
157 e_z[1] = eps_z_i(1, 1);
158 e_z[2] = eps_z_i(2, 2);
159 e_z[3] = 2 * eps_z_i(1, 2);
160 e_z[4] = 2 * eps_z_i(0, 2);
161 e_z[5] = 2 * eps_z_i(0, 1);
162
163 Eigen::Matrix3d sigma_x;
167
168 Eigen::Matrix3d sigma_y;
172
173 Eigen::Matrix3d sigma_z;
177
178 res_k(0) = (sigma_x * eps_x_j).trace();
179 res_k(3) = (sigma_x * eps_y_j).trace();
180 res_k(6) = (sigma_x * eps_z_j).trace();
181
182 res_k(1) = (sigma_y * eps_x_j).trace();
183 res_k(4) = (sigma_y * eps_y_j).trace();
184 res_k(7) = (sigma_y * eps_z_j).trace();
185
186 res_k(2) = (sigma_z * eps_x_j).trace();
187 res_k(5) = (sigma_z * eps_y_j).trace();
188 res_k(8) = (sigma_z * eps_z_j).trace();
189 }
190
191 res += res_k * data.da(k);
192 }
193
194 return res;
195 }
196
198 const OutputData &data,
199 const int all_size,
200 const ElasticityTensorType &type,
201 Eigen::MatrixXd &all,
202 const std::function<Eigen::MatrixXd(const Eigen::MatrixXd &)> &fun) const
203 {
204 const auto &displacement = data.fun;
205 const auto &local_pts = data.local_pts;
206 const auto &bs = data.bs;
207 const auto &gbs = data.gbs;
208 const auto el_id = data.el_id;
209
210 Eigen::MatrixXd displacement_grad(size(), size());
211
212 assert(displacement.cols() == 1);
213
214 all.resize(local_pts.rows(), all_size);
215
217 vals.compute(el_id, size() == 3, local_pts, bs, gbs);
218
219 for (long p = 0; p < local_pts.rows(); ++p)
220 {
221 compute_diplacement_grad(size(), bs, vals, local_pts, p, displacement, displacement_grad);
222
223 if (type == ElasticityTensorType::F)
224 {
225 all.row(p) = fun(displacement_grad + Eigen::MatrixXd::Identity(size(), size()));
226 continue;
227 }
228
229 Eigen::MatrixXd strain = (displacement_grad + displacement_grad.transpose()) / 2;
230 Eigen::MatrixXd sigma(size(), size());
231
232 if (size() == 2)
233 {
234 std::array<double, 3> eps;
235 eps[0] = strain(0, 0);
236 eps[1] = strain(1, 1);
237 eps[2] = 2 * strain(0, 1);
238
241 }
242 else
243 {
244 std::array<double, 6> eps;
245 eps[0] = strain(0, 0);
246 eps[1] = strain(1, 1);
247 eps[2] = strain(2, 2);
248 eps[3] = 2 * strain(1, 2);
249 eps[4] = 2 * strain(0, 2);
250 eps[5] = 2 * strain(0, 1);
251
255 }
256
257 if (type == ElasticityTensorType::PK1)
258 sigma = pk1_from_cauchy(sigma, displacement_grad + Eigen::MatrixXd::Identity(size(), size()));
259 else if (type == ElasticityTensorType::PK2)
260 sigma = pk2_from_cauchy(sigma, displacement_grad + Eigen::MatrixXd::Identity(size(), size()));
261
262 all.row(p) = fun(sigma);
263 }
264 }
265
266 Eigen::Matrix<double, Eigen::Dynamic, 1, 0, 3, 1>
268 {
269 assert(pt.size() == size());
270 Eigen::Matrix<double, Eigen::Dynamic, 1, 0, 3, 1> res;
271
272 if (size() == 2)
274 else if (size() == 3)
276 else
277 assert(false);
278
279 return res;
280 }
281
282 std::map<std::string, Assembler::ParamFunc> HookeLinearElasticity::parameters() const
283 {
284 std::map<std::string, ParamFunc> res;
285 const auto &elast_tensor = elasticity_tensor();
286 const int size = this->size() == 2 ? 3 : 6;
287
288 for (int i = 0; i < size; ++i)
289 {
290 for (int j = i; j < size; ++j)
291 {
292 res[fmt::format("C_{}{}", i, j)] = [&elast_tensor, i, j](const RowVectorNd &, const RowVectorNd &, double, int) {
293 return elast_tensor(i, j);
294 };
295 }
296 }
297 return res;
298 }
299
301 {
302 return compute_energy_aux<double>(data);
303 }
304
306 {
307 const int n_bases = data.vals.basis_values.size();
309 size(), n_bases, data,
310 [&](const NonLinearAssemblerData &data) { return compute_energy_aux<DScalar1<double, Eigen::Matrix<double, 6, 1>>>(data); },
311 [&](const NonLinearAssemblerData &data) { return compute_energy_aux<DScalar1<double, Eigen::Matrix<double, 8, 1>>>(data); },
312 [&](const NonLinearAssemblerData &data) { return compute_energy_aux<DScalar1<double, Eigen::Matrix<double, 12, 1>>>(data); },
313 [&](const NonLinearAssemblerData &data) { return compute_energy_aux<DScalar1<double, Eigen::Matrix<double, 18, 1>>>(data); },
314 [&](const NonLinearAssemblerData &data) { return compute_energy_aux<DScalar1<double, Eigen::Matrix<double, 24, 1>>>(data); },
315 [&](const NonLinearAssemblerData &data) { return compute_energy_aux<DScalar1<double, Eigen::Matrix<double, 30, 1>>>(data); },
316 [&](const NonLinearAssemblerData &data) { return compute_energy_aux<DScalar1<double, Eigen::Matrix<double, 60, 1>>>(data); },
317 [&](const NonLinearAssemblerData &data) { return compute_energy_aux<DScalar1<double, Eigen::Matrix<double, 81, 1>>>(data); },
318 [&](const NonLinearAssemblerData &data) { return compute_energy_aux<DScalar1<double, Eigen::Matrix<double, Eigen::Dynamic, 1, 0, SMALL_N, 1>>>(data); },
319 [&](const NonLinearAssemblerData &data) { return compute_energy_aux<DScalar1<double, Eigen::Matrix<double, Eigen::Dynamic, 1, 0, BIG_N, 1>>>(data); },
320 [&](const NonLinearAssemblerData &data) { return compute_energy_aux<DScalar1<double, Eigen::VectorXd>>(data); });
321 }
322
324 {
325 const int n_bases = data.vals.basis_values.size();
327 size(), n_bases, data,
328 [&](const NonLinearAssemblerData &data) { return compute_energy_aux<DScalar2<double, Eigen::Matrix<double, 6, 1>, Eigen::Matrix<double, 6, 6>>>(data); },
329 [&](const NonLinearAssemblerData &data) { return compute_energy_aux<DScalar2<double, Eigen::Matrix<double, 8, 1>, Eigen::Matrix<double, 8, 8>>>(data); },
330 [&](const NonLinearAssemblerData &data) { return compute_energy_aux<DScalar2<double, Eigen::Matrix<double, 12, 1>, Eigen::Matrix<double, 12, 12>>>(data); },
331 [&](const NonLinearAssemblerData &data) { return compute_energy_aux<DScalar2<double, Eigen::Matrix<double, 18, 1>, Eigen::Matrix<double, 18, 18>>>(data); },
332 [&](const NonLinearAssemblerData &data) { return compute_energy_aux<DScalar2<double, Eigen::Matrix<double, 24, 1>, Eigen::Matrix<double, 24, 24>>>(data); },
333 [&](const NonLinearAssemblerData &data) { return compute_energy_aux<DScalar2<double, Eigen::Matrix<double, 30, 1>, Eigen::Matrix<double, 30, 30>>>(data); },
334 [&](const NonLinearAssemblerData &data) { return compute_energy_aux<DScalar2<double, Eigen::Matrix<double, 60, 1>, Eigen::Matrix<double, 60, 60>>>(data); },
335 [&](const NonLinearAssemblerData &data) { return compute_energy_aux<DScalar2<double, Eigen::Matrix<double, 81, 1>, Eigen::Matrix<double, 81, 81>>>(data); },
336 [&](const NonLinearAssemblerData &data) { return compute_energy_aux<DScalar2<double, Eigen::Matrix<double, Eigen::Dynamic, 1, 0, SMALL_N, 1>, Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, 0, SMALL_N, SMALL_N>>>(data); },
337 [&](const NonLinearAssemblerData &data) { return compute_energy_aux<DScalar2<double, Eigen::VectorXd, Eigen::MatrixXd>>(data); });
338 }
339
340 template <typename T>
342 {
343 typedef Eigen::Matrix<T, Eigen::Dynamic, 1> AutoDiffVect;
344 typedef Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic, 0, 3, 3> AutoDiffGradMat;
345
346 AutoDiffVect local_disp;
347 get_local_disp(data, size(), local_disp);
348
349 AutoDiffGradMat disp_grad(size(), size());
350
351 T energy = T(0.0);
352
353 const int n_pts = data.da.size();
354 for (long p = 0; p < n_pts; ++p)
355 {
356 compute_disp_grad_at_quad(data, local_disp, p, size(), disp_grad);
357
358 AutoDiffGradMat strain = strain_from_disp_grad(disp_grad);
359 AutoDiffGradMat stress_tensor(size(), size());
360
361 if (size() == 2)
362 {
363 std::array<T, 3> eps;
364 eps[0] = strain(0, 0);
365 eps[1] = strain(1, 1);
366 eps[2] = 2 * strain(0, 1);
367
368 stress_tensor << stress(elasticity_tensor_, eps, 0), stress(elasticity_tensor_, eps, 2),
369 stress(elasticity_tensor_, eps, 2), stress(elasticity_tensor_, eps, 1);
370 }
371 else
372 {
373 std::array<T, 6> eps;
374 eps[0] = strain(0, 0);
375 eps[1] = strain(1, 1);
376 eps[2] = strain(2, 2);
377 eps[3] = 2 * strain(1, 2);
378 eps[4] = 2 * strain(0, 2);
379 eps[5] = 2 * strain(0, 1);
380
381 stress_tensor << stress(elasticity_tensor_, eps, 0), stress(elasticity_tensor_, eps, 5), stress(elasticity_tensor_, eps, 4),
382 stress(elasticity_tensor_, eps, 5), stress(elasticity_tensor_, eps, 1), stress(elasticity_tensor_, eps, 3),
383 stress(elasticity_tensor_, eps, 4), stress(elasticity_tensor_, eps, 3), stress(elasticity_tensor_, eps, 2);
384 }
385
386 energy += (stress_tensor * strain).trace() * data.da(p);
387 }
388
389 return energy * 0.5;
390 }
391} // namespace polyfem::assembler
ElementAssemblyValues vals
Definition Assembler.cpp:22
std::vector< Eigen::Triplet< double > > entries
std::string stress() const
Definition Units.hpp:27
virtual void set_size(const int size)
Definition Assembler.hpp:64
void set_from_entries(const std::vector< double > &entries, const std::string &stress_unit)
void set_from_young_poisson(const double young, const double poisson, const std::string &stress_unit)
double compute_stress(const std::array< double, DIM > &strain, const int j) const
void set_from_lambda_mu(const double lambda, const double mu, const std::string &stress_unit)
stores per element basis values at given quadrature points and geometric mapping
void compute(const int el_index, const bool is_volume, const Eigen::MatrixXd &pts, const basis::ElementBases &basis, const basis::ElementBases &gbasis)
computes the per element values at the local (ref el) points (pts) sets basis_values,...
std::vector< Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic, 0, 3, 3 > > jac_it
Eigen::Matrix< double, Eigen::Dynamic, 1, 0, 9, 1 > assemble(const LinearAssemblerData &data) const override
local assembly function that defines the bilinear form (LHS) computes and returns a single local stif...
double compute_energy(const NonLinearAssemblerData &data) const override
std::map< std::string, ParamFunc > parameters() const override
const ElasticityTensor & elasticity_tensor() const
void add_multimaterial(const int index, const json &params, const Units &units) override
VectorNd compute_rhs(const AutodiffHessianPt &pt) const override
Eigen::MatrixXd assemble_hessian(const NonLinearAssemblerData &data) const override
T compute_energy_aux(const NonLinearAssemblerData &data) const
Eigen::VectorXd assemble_gradient(const NonLinearAssemblerData &data) const override
void assign_stress_tensor(const OutputData &data, const int all_size, const ElasticityTensorType &type, Eigen::MatrixXd &all, const std::function< Eigen::MatrixXd(const Eigen::MatrixXd &)> &fun) const override
const ElementAssemblyValues & vals
stores the evaluation for that element
const QuadratureVector & da
contains both the quadrature weight and the change of metric in the integral
const basis::ElementBases & bs
const Eigen::MatrixXd & fun
const Eigen::MatrixXd & local_pts
const basis::ElementBases & gbs
Used for test only.
Eigen::Matrix< T, Eigen::Dynamic, Eigen::Dynamic, 0, 3, 3 > AutoDiffGradMat
Eigen::Matrix< T, Eigen::Dynamic, 1 > AutoDiffVect
void hooke_3d_function(const AutodiffHessianPt &pt, const assembler::ElasticityTensor &C, Eigen::Matrix< double, Eigen::Dynamic, 1, 0, 3, 1 > &res)
void hooke_2d_function(const AutodiffHessianPt &pt, const assembler::ElasticityTensor &C, Eigen::Matrix< double, Eigen::Dynamic, 1, 0, 3, 1 > &res)
Eigen::MatrixXd pk2_from_cauchy(const Eigen::MatrixXd &stress, const Eigen::MatrixXd &F)
void compute_disp_grad_at_quad(const assembler::NonLinearAssemblerData &data, const AutoDiffVect &local_disp, const int p, const int size, AutoDiffGradMat &def_grad)
Eigen::Matrix< AutodiffScalarHessian, Eigen::Dynamic, 1, 0, 3, 1 > AutodiffHessianPt
Eigen::MatrixXd hessian_from_energy(const int size, const int n_bases, const assembler::NonLinearAssemblerData &data, const std::function< DScalar2< double, Eigen::Matrix< double, 6, 1 >, Eigen::Matrix< double, 6, 6 > >(const assembler::NonLinearAssemblerData &)> &fun6, const std::function< DScalar2< double, Eigen::Matrix< double, 8, 1 >, Eigen::Matrix< double, 8, 8 > >(const assembler::NonLinearAssemblerData &)> &fun8, const std::function< DScalar2< double, Eigen::Matrix< double, 12, 1 >, Eigen::Matrix< double, 12, 12 > >(const assembler::NonLinearAssemblerData &)> &fun12, const std::function< DScalar2< double, Eigen::Matrix< double, 18, 1 >, Eigen::Matrix< double, 18, 18 > >(const assembler::NonLinearAssemblerData &)> &fun18, const std::function< DScalar2< double, Eigen::Matrix< double, 24, 1 >, Eigen::Matrix< double, 24, 24 > >(const assembler::NonLinearAssemblerData &)> &fun24, const std::function< DScalar2< double, Eigen::Matrix< double, 30, 1 >, Eigen::Matrix< double, 30, 30 > >(const assembler::NonLinearAssemblerData &)> &fun30, const std::function< DScalar2< double, Eigen::Matrix< double, 60, 1 >, Eigen::Matrix< double, 60, 60 > >(const assembler::NonLinearAssemblerData &)> &fun60, const std::function< DScalar2< double, Eigen::Matrix< double, 81, 1 >, Eigen::Matrix< double, 81, 81 > >(const assembler::NonLinearAssemblerData &)> &fun81, const std::function< DScalar2< double, Eigen::Matrix< double, Eigen::Dynamic, 1, 0, SMALL_N, 1 >, Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic, 0, SMALL_N, SMALL_N > >(const assembler::NonLinearAssemblerData &)> &funN, const std::function< DScalar2< double, Eigen::VectorXd, Eigen::MatrixXd >(const assembler::NonLinearAssemblerData &)> &funn)
nlohmann::json json
Definition Common.hpp:9
void compute_diplacement_grad(const int size, const ElementAssemblyValues &vals, const Eigen::MatrixXd &local_pts, const int p, const Eigen::MatrixXd &displacement, Eigen::MatrixXd &displacement_grad)
void get_local_disp(const assembler::NonLinearAssemblerData &data, const int size, AutoDiffVect &local_disp)
Eigen::MatrixXd pk1_from_cauchy(const Eigen::MatrixXd &stress, const Eigen::MatrixXd &F)
Eigen::Matrix< double, 1, Eigen::Dynamic, Eigen::RowMajor, 1, 3 > RowVectorNd
Definition Types.hpp:13
Eigen::VectorXd gradient_from_energy(const int size, const int n_bases, const assembler::NonLinearAssemblerData &data, const std::function< DScalar1< double, Eigen::Matrix< double, 6, 1 > >(const assembler::NonLinearAssemblerData &)> &fun6, const std::function< DScalar1< double, Eigen::Matrix< double, 8, 1 > >(const assembler::NonLinearAssemblerData &)> &fun8, const std::function< DScalar1< double, Eigen::Matrix< double, 12, 1 > >(const assembler::NonLinearAssemblerData &)> &fun12, const std::function< DScalar1< double, Eigen::Matrix< double, 18, 1 > >(const assembler::NonLinearAssemblerData &)> &fun18, const std::function< DScalar1< double, Eigen::Matrix< double, 24, 1 > >(const assembler::NonLinearAssemblerData &)> &fun24, const std::function< DScalar1< double, Eigen::Matrix< double, 30, 1 > >(const assembler::NonLinearAssemblerData &)> &fun30, const std::function< DScalar1< double, Eigen::Matrix< double, 60, 1 > >(const assembler::NonLinearAssemblerData &)> &fun60, const std::function< DScalar1< double, Eigen::Matrix< double, 81, 1 > >(const assembler::NonLinearAssemblerData &)> &fun81, const std::function< DScalar1< double, Eigen::Matrix< double, Eigen::Dynamic, 1, 0, SMALL_N, 1 > >(const assembler::NonLinearAssemblerData &)> &funN, const std::function< DScalar1< double, Eigen::Matrix< double, Eigen::Dynamic, 1, 0, BIG_N, 1 > >(const assembler::NonLinearAssemblerData &)> &funBigN, const std::function< DScalar1< double, Eigen::VectorXd >(const assembler::NonLinearAssemblerData &)> &funn)