51 Assert(values.size() == this->n() || values.empty(),
53 Assert(grads.size() == this->n() || grads.empty(),
55 Assert(grad_grads.size() == this->n() || grad_grads.empty(),
57 Assert(third_derivatives.size() == this->n() || third_derivatives.empty(),
59 Assert(fourth_derivatives.size() == this->n() || fourth_derivatives.empty(),
63 (void)third_derivatives;
65 (void)fourth_derivatives;
68 const unsigned int n_sub = raviart_thomas_space.n();
69 const unsigned int my_degree = this->degree();
75 static std::mutex mutex;
76 std::lock_guard<std::mutex> lock(mutex);
78 static std::vector<Tensor<1, dim>> p_values;
79 static std::vector<Tensor<2, dim>> p_grads;
80 static std::vector<Tensor<3, dim>> p_grad_grads;
81 static std::vector<Tensor<4, dim>> p_third_derivatives;
82 static std::vector<Tensor<5, dim>> p_fourth_derivatives;
84 p_values.resize((values.empty()) ? 0 : n_sub);
85 p_grads.resize((grads.empty()) ? 0 : n_sub);
86 p_grad_grads.resize((grad_grads.empty()) ? 0 : n_sub);
89 raviart_thomas_space.evaluate(unit_point,
94 p_fourth_derivatives);
95 for (
unsigned int i = 0; i < p_values.size(); ++i)
96 values[i] = p_values[i];
97 for (
unsigned int i = 0; i < p_grads.size(); ++i)
98 grads[i] = p_grads[i];
99 for (
unsigned int i = 0; i < p_grad_grads.size(); ++i)
100 grad_grads[i] = p_grad_grads[i];
105 const unsigned int n_derivatives = 3;
106 double monoval_plus[dim][n_derivatives + 1];
107 double monoval_i[dim][n_derivatives + 1];
110 if constexpr (dim <= 1)
116 unsigned int start = n_sub;
117 if constexpr (dim == 2)
127 for (
unsigned int d = 0; d < dim; ++d)
128 monomials[my_degree + 1].value(unit_point[d],
132 for (
unsigned int i = 0; i <= my_degree; ++i, ++start)
134 for (
unsigned int d = 0; d < dim; ++d)
135 monomials[i].value(unit_point[d], n_derivatives, monoval_i[d]);
137 if (values.size() != 0)
139 values[start][0] = monoval_i[0][0] * monoval_plus[1][1];
140 values[start][1] = -monoval_i[0][1] * monoval_plus[1][0];
142 values[start + my_degree + 1][0] =
143 -monoval_plus[0][0] * monoval_i[1][1];
144 values[start + my_degree + 1][1] =
145 monoval_plus[0][1] * monoval_i[1][0];
148 if (grads.size() != 0)
150 grads[start][0][0] = monoval_i[0][1] * monoval_plus[1][1];
151 grads[start][0][1] = monoval_i[0][0] * monoval_plus[1][2];
152 grads[start][1][0] = -monoval_i[0][2] * monoval_plus[1][0];
153 grads[start][1][1] = -monoval_i[0][1] * monoval_plus[1][1];
155 grads[start + my_degree + 1][0][0] =
156 -monoval_plus[0][1] * monoval_i[1][1];
157 grads[start + my_degree + 1][0][1] =
158 -monoval_plus[0][0] * monoval_i[1][2];
159 grads[start + my_degree + 1][1][0] =
160 monoval_plus[0][2] * monoval_i[1][0];
161 grads[start + my_degree + 1][1][1] =
162 monoval_plus[0][1] * monoval_i[1][1];
165 if (grad_grads.size() != 0)
167 grad_grads[start][0][0][0] = monoval_i[0][2] * monoval_plus[1][1];
168 grad_grads[start][0][0][1] = monoval_i[0][1] * monoval_plus[1][2];
169 grad_grads[start][0][1][0] = monoval_i[0][1] * monoval_plus[1][2];
170 grad_grads[start][0][1][1] = monoval_i[0][0] * monoval_plus[1][3];
171 grad_grads[start][1][0][0] =
172 -monoval_i[0][3] * monoval_plus[1][0];
173 grad_grads[start][1][0][1] =
174 -monoval_i[0][2] * monoval_plus[1][1];
175 grad_grads[start][1][1][0] =
176 -monoval_i[0][2] * monoval_plus[1][1];
177 grad_grads[start][1][1][1] =
178 -monoval_i[0][1] * monoval_plus[1][2];
180 grad_grads[start + my_degree + 1][0][0][0] =
181 -monoval_plus[0][2] * monoval_i[1][1];
182 grad_grads[start + my_degree + 1][0][0][1] =
183 -monoval_plus[0][1] * monoval_i[1][2];
184 grad_grads[start + my_degree + 1][0][1][0] =
185 -monoval_plus[0][1] * monoval_i[1][2];
186 grad_grads[start + my_degree + 1][0][1][1] =
187 -monoval_plus[0][0] * monoval_i[1][3];
188 grad_grads[start + my_degree + 1][1][0][0] =
189 monoval_plus[0][3] * monoval_i[1][0];
190 grad_grads[start + my_degree + 1][1][0][1] =
191 monoval_plus[0][2] * monoval_i[1][1];
192 grad_grads[start + my_degree + 1][1][1][0] =
193 monoval_plus[0][2] * monoval_i[1][1];
194 grad_grads[start + my_degree + 1][1][1][1] =
195 monoval_plus[0][1] * monoval_i[1][2];
200 else if constexpr (dim == 3)
202 double monoval[dim][n_derivatives + 1];
203 double monoval_j[dim][n_derivatives + 1];
204 double monoval_jplus[dim][n_derivatives + 1];
216 for (
unsigned int d = 0; d < dim; ++d)
218 monomials[my_degree + 1].value(unit_point[d],
221 monomials[my_degree].value(unit_point[d], n_derivatives, monoval[d]);
224 const unsigned int n_curls = (my_degree + 1) * (2 * my_degree + 1);
226 for (
unsigned int i = 0; i <= my_degree; ++i)
228 for (
unsigned int d = 0; d < dim; ++d)
229 monomials[i].value(unit_point[d], n_derivatives, monoval_i[d]);
231 for (
unsigned int j = 0; j <= my_degree; ++j)
233 for (
unsigned int d = 0; d < dim; ++d)
235 monomials[j].value(unit_point[d],
238 monomials[j + 1].value(unit_point[d],
243 if (values.size() != 0)
245 values[start][0] = monoval_i[0][0] * monoval_j[1][0] *
247 static_cast<double>(j + my_degree + 2);
249 -monoval_i[0][1] * monoval_jplus[1][0] * monoval[2][0];
251 -monoval_i[0][1] * monoval_j[1][0] * monoval_plus[2][0];
253 values[start + n_curls][0] =
254 -monoval_jplus[0][0] * monoval_i[1][1] * monoval[2][0];
255 values[start + n_curls][1] =
256 monoval_j[0][0] * monoval_i[1][0] * monoval[2][0] *
257 static_cast<double>(j + my_degree + 2);
258 values[start + n_curls][2] =
259 -monoval_j[0][0] * monoval_i[1][1] * monoval_plus[2][0];
261 values[start + 2 * n_curls][0] =
262 -monoval_jplus[0][0] * monoval[1][0] * monoval_i[2][1];
263 values[start + 2 * n_curls][1] =
264 -monoval_j[0][0] * monoval_plus[1][0] * monoval_i[2][1];
265 values[start + 2 * n_curls][2] =
266 monoval_j[0][0] * monoval[1][0] * monoval_i[2][0] *
267 static_cast<double>(j + my_degree + 2);
273 values[start + 1][0] =
274 monoval_i[0][0] * monoval[1][0] * monoval_j[2][0] *
275 static_cast<double>(j + my_degree + 2);
276 values[start + 1][1] =
277 -monoval_i[0][1] * monoval_plus[1][0] * monoval_j[2][0];
278 values[start + 1][2] =
279 -monoval_i[0][1] * monoval[1][0] * monoval_jplus[2][0];
281 values[start + n_curls + 1][0] =
282 -monoval_plus[0][0] * monoval_i[1][1] * monoval_j[2][0];
283 values[start + n_curls + 1][1] =
284 monoval[0][0] * monoval_i[1][0] * monoval_j[2][0] *
285 static_cast<double>(j + my_degree + 2);
286 values[start + n_curls + 1][2] =
287 -monoval[0][0] * monoval_i[1][1] * monoval_jplus[2][0];
289 values[start + 2 * n_curls + 1][0] =
290 -monoval_plus[0][0] * monoval_j[1][0] * monoval_i[2][1];
291 values[start + 2 * n_curls + 1][1] =
292 -monoval[0][0] * monoval_jplus[1][0] * monoval_i[2][1];
293 values[start + 2 * n_curls + 1][2] =
294 monoval[0][0] * monoval_j[1][0] * monoval_i[2][0] *
295 static_cast<double>(j + my_degree + 2);
299 if (grads.size() != 0)
301 grads[start][0][0] = monoval_i[0][1] * monoval_j[1][0] *
303 static_cast<double>(j + my_degree + 2);
304 grads[start][0][1] = monoval_i[0][0] * monoval_j[1][1] *
306 static_cast<double>(j + my_degree + 2);
307 grads[start][0][2] = monoval_i[0][0] * monoval_j[1][0] *
309 static_cast<double>(j + my_degree + 2);
311 -monoval_i[0][2] * monoval_jplus[1][0] * monoval[2][0];
313 -monoval_i[0][1] * monoval_jplus[1][1] * monoval[2][0];
315 -monoval_i[0][1] * monoval_jplus[1][0] * monoval[2][1];
317 -monoval_i[0][2] * monoval_j[1][0] * monoval_plus[2][0];
319 -monoval_i[0][1] * monoval_j[1][1] * monoval_plus[2][0];
321 -monoval_i[0][1] * monoval_j[1][0] * monoval_plus[2][1];
323 grads[start + n_curls][0][0] =
324 -monoval_jplus[0][1] * monoval_i[1][1] * monoval[2][0];
325 grads[start + n_curls][0][1] =
326 -monoval_jplus[0][0] * monoval_i[1][2] * monoval[2][0];
327 grads[start + n_curls][0][2] =
328 -monoval_jplus[0][0] * monoval_i[1][1] * monoval[2][1];
329 grads[start + n_curls][1][0] =
330 monoval_j[0][1] * monoval_i[1][0] * monoval[2][0] *
331 static_cast<double>(j + my_degree + 2);
332 grads[start + n_curls][1][1] =
333 monoval_j[0][0] * monoval_i[1][1] * monoval[2][0] *
334 static_cast<double>(j + my_degree + 2);
335 grads[start + n_curls][1][2] =
336 monoval_j[0][0] * monoval_i[1][0] * monoval[2][1] *
337 static_cast<double>(j + my_degree + 2);
338 grads[start + n_curls][2][0] =
339 -monoval_j[0][1] * monoval_i[1][1] * monoval_plus[2][0];
340 grads[start + n_curls][2][1] =
341 -monoval_j[0][0] * monoval_i[1][2] * monoval_plus[2][0];
342 grads[start + n_curls][2][2] =
343 -monoval_j[0][0] * monoval_i[1][1] * monoval_plus[2][1];
345 grads[start + 2 * n_curls][0][0] =
346 -monoval_jplus[0][1] * monoval[1][0] * monoval_i[2][1];
347 grads[start + 2 * n_curls][0][1] =
348 -monoval_jplus[0][0] * monoval[1][1] * monoval_i[2][1];
349 grads[start + 2 * n_curls][0][2] =
350 -monoval_jplus[0][0] * monoval[1][0] * monoval_i[2][2];
351 grads[start + 2 * n_curls][1][0] =
352 -monoval_j[0][1] * monoval_plus[1][0] * monoval_i[2][1];
353 grads[start + 2 * n_curls][1][1] =
354 -monoval_j[0][0] * monoval_plus[1][1] * monoval_i[2][1];
355 grads[start + 2 * n_curls][1][2] =
356 -monoval_j[0][0] * monoval_plus[1][0] * monoval_i[2][2];
357 grads[start + 2 * n_curls][2][0] =
358 monoval_j[0][1] * monoval[1][0] * monoval_i[2][0] *
359 static_cast<double>(j + my_degree + 2);
360 grads[start + 2 * n_curls][2][1] =
361 monoval_j[0][0] * monoval[1][1] * monoval_i[2][0] *
362 static_cast<double>(j + my_degree + 2);
363 grads[start + 2 * n_curls][2][2] =
364 monoval_j[0][0] * monoval[1][0] * monoval_i[2][1] *
365 static_cast<double>(j + my_degree + 2);
369 grads[start + 1][0][0] =
370 monoval_i[0][1] * monoval[1][0] * monoval_j[2][0] *
371 static_cast<double>(j + my_degree + 2);
372 grads[start + 1][0][1] =
373 monoval_i[0][0] * monoval[1][1] * monoval_j[2][0] *
374 static_cast<double>(j + my_degree + 2);
375 grads[start + 1][0][2] =
376 monoval_i[0][0] * monoval[1][0] * monoval_j[2][1] *
377 static_cast<double>(j + my_degree + 2);
378 grads[start + 1][1][0] =
379 -monoval_i[0][2] * monoval_plus[1][0] * monoval_j[2][0];
380 grads[start + 1][1][1] =
381 -monoval_i[0][1] * monoval_plus[1][1] * monoval_j[2][0];
382 grads[start + 1][1][2] =
383 -monoval_i[0][1] * monoval_plus[1][0] * monoval_j[2][1];
384 grads[start + 1][2][0] =
385 -monoval_i[0][2] * monoval[1][0] * monoval_jplus[2][0];
386 grads[start + 1][2][1] =
387 -monoval_i[0][1] * monoval[1][1] * monoval_jplus[2][0];
388 grads[start + 1][2][2] =
389 -monoval_i[0][1] * monoval[1][0] * monoval_jplus[2][1];
391 grads[start + n_curls + 1][0][0] =
392 -monoval_plus[0][1] * monoval_i[1][1] * monoval_j[2][0];
393 grads[start + n_curls + 1][0][1] =
394 -monoval_plus[0][0] * monoval_i[1][2] * monoval_j[2][0];
395 grads[start + n_curls + 1][0][2] =
396 -monoval_plus[0][0] * monoval_i[1][1] * monoval_j[2][1];
397 grads[start + n_curls + 1][1][0] =
398 monoval[0][1] * monoval_i[1][0] * monoval_j[2][0] *
399 static_cast<double>(j + my_degree + 2);
400 grads[start + n_curls + 1][1][1] =
401 monoval[0][0] * monoval_i[1][1] * monoval_j[2][0] *
402 static_cast<double>(j + my_degree + 2);
403 grads[start + n_curls + 1][1][2] =
404 monoval[0][0] * monoval_i[1][0] * monoval_j[2][1] *
405 static_cast<double>(j + my_degree + 2);
406 grads[start + n_curls + 1][2][0] =
407 -monoval[0][1] * monoval_i[1][1] * monoval_jplus[2][0];
408 grads[start + n_curls + 1][2][1] =
409 -monoval[0][0] * monoval_i[1][2] * monoval_jplus[2][0];
410 grads[start + n_curls + 1][2][2] =
411 -monoval[0][0] * monoval_i[1][1] * monoval_jplus[2][1];
413 grads[start + 2 * n_curls + 1][0][0] =
414 -monoval_plus[0][1] * monoval_j[1][0] * monoval_i[2][1];
415 grads[start + 2 * n_curls + 1][0][1] =
416 -monoval_plus[0][0] * monoval_j[1][1] * monoval_i[2][1];
417 grads[start + 2 * n_curls + 1][0][2] =
418 -monoval_plus[0][0] * monoval_j[1][0] * monoval_i[2][2];
419 grads[start + 2 * n_curls + 1][1][0] =
420 -monoval[0][1] * monoval_jplus[1][0] * monoval_i[2][1];
421 grads[start + 2 * n_curls + 1][1][1] =
422 -monoval[0][0] * monoval_jplus[1][1] * monoval_i[2][1];
423 grads[start + 2 * n_curls + 1][1][2] =
424 -monoval[0][0] * monoval_jplus[1][0] * monoval_i[2][2];
425 grads[start + 2 * n_curls + 1][2][0] =
426 monoval[0][1] * monoval_j[1][0] * monoval_i[2][0] *
427 static_cast<double>(j + my_degree + 2);
428 grads[start + 2 * n_curls + 1][2][1] =
429 monoval[0][0] * monoval_j[1][1] * monoval_i[2][0] *
430 static_cast<double>(j + my_degree + 2);
431 grads[start + 2 * n_curls + 1][2][2] =
432 monoval[0][0] * monoval_j[1][0] * monoval_i[2][1] *
433 static_cast<double>(j + my_degree + 2);
437 if (grad_grads.size() != 0)
439 grad_grads[start][0][0][0] =
440 monoval_i[0][2] * monoval_j[1][0] * monoval[2][0] *
441 static_cast<double>(j + my_degree + 2);
442 grad_grads[start][0][0][1] =
443 monoval_i[0][1] * monoval_j[1][1] * monoval[2][0] *
444 static_cast<double>(j + my_degree + 2);
445 grad_grads[start][0][0][2] =
446 monoval_i[0][1] * monoval_j[1][0] * monoval[2][1] *
447 static_cast<double>(j + my_degree + 2);
448 grad_grads[start][0][1][0] =
449 monoval_i[0][1] * monoval_j[1][1] * monoval[2][0] *
450 static_cast<double>(j + my_degree + 2);
451 grad_grads[start][0][1][1] =
452 monoval_i[0][0] * monoval_j[1][2] * monoval[2][0] *
453 static_cast<double>(j + my_degree + 2);
454 grad_grads[start][0][1][2] =
455 monoval_i[0][0] * monoval_j[1][1] * monoval[2][1] *
456 static_cast<double>(j + my_degree + 2);
457 grad_grads[start][0][2][0] =
458 monoval_i[0][1] * monoval_j[1][0] * monoval[2][1] *
459 static_cast<double>(j + my_degree + 2);
460 grad_grads[start][0][2][1] =
461 monoval_i[0][0] * monoval_j[1][1] * monoval[2][1] *
462 static_cast<double>(j + my_degree + 2);
463 grad_grads[start][0][2][2] =
464 monoval_i[0][0] * monoval_j[1][0] * monoval[2][2] *
465 static_cast<double>(j + my_degree + 2);
466 grad_grads[start][1][0][0] =
467 -monoval_i[0][3] * monoval_jplus[1][0] * monoval[2][0];
468 grad_grads[start][1][0][1] =
469 -monoval_i[0][2] * monoval_jplus[1][1] * monoval[2][0];
470 grad_grads[start][1][0][2] =
471 -monoval_i[0][2] * monoval_jplus[1][0] * monoval[2][1];
472 grad_grads[start][1][1][0] =
473 -monoval_i[0][2] * monoval_jplus[1][1] * monoval[2][0];
474 grad_grads[start][1][1][1] =
475 -monoval_i[0][1] * monoval_jplus[1][2] * monoval[2][0];
476 grad_grads[start][1][1][2] =
477 -monoval_i[0][1] * monoval_jplus[1][1] * monoval[2][1];
478 grad_grads[start][1][2][0] =
479 -monoval_i[0][2] * monoval_jplus[1][0] * monoval[2][1];
480 grad_grads[start][1][2][1] =
481 -monoval_i[0][1] * monoval_jplus[1][1] * monoval[2][1];
482 grad_grads[start][1][2][2] =
483 -monoval_i[0][1] * monoval_jplus[1][0] * monoval[2][2];
484 grad_grads[start][2][0][0] =
485 -monoval_i[0][3] * monoval_j[1][0] * monoval_plus[2][0];
486 grad_grads[start][2][0][1] =
487 -monoval_i[0][2] * monoval_j[1][1] * monoval_plus[2][0];
488 grad_grads[start][2][0][2] =
489 -monoval_i[0][2] * monoval_j[1][0] * monoval_plus[2][1];
490 grad_grads[start][2][1][0] =
491 -monoval_i[0][2] * monoval_j[1][1] * monoval_plus[2][0];
492 grad_grads[start][2][1][1] =
493 -monoval_i[0][1] * monoval_j[1][2] * monoval_plus[2][0];
494 grad_grads[start][2][1][2] =
495 -monoval_i[0][1] * monoval_j[1][1] * monoval_plus[2][1];
496 grad_grads[start][2][2][0] =
497 -monoval_i[0][2] * monoval_j[1][0] * monoval_plus[2][1];
498 grad_grads[start][2][2][1] =
499 -monoval_i[0][1] * monoval_j[1][1] * monoval_plus[2][1];
500 grad_grads[start][2][2][2] =
501 -monoval_i[0][1] * monoval_j[1][0] * monoval_plus[2][2];
503 grad_grads[start + n_curls][0][0][0] =
504 -monoval_jplus[0][2] * monoval_i[1][1] * monoval[2][0];
505 grad_grads[start + n_curls][0][0][1] =
506 -monoval_jplus[0][1] * monoval_i[1][2] * monoval[2][0];
507 grad_grads[start + n_curls][0][0][2] =
508 -monoval_jplus[0][1] * monoval_i[1][1] * monoval[2][1];
509 grad_grads[start + n_curls][0][1][0] =
510 -monoval_jplus[0][1] * monoval_i[1][2] * monoval[2][0];
511 grad_grads[start + n_curls][0][1][1] =
512 -monoval_jplus[0][0] * monoval_i[1][3] * monoval[2][0];
513 grad_grads[start + n_curls][0][1][2] =
514 -monoval_jplus[0][0] * monoval_i[1][2] * monoval[2][1];
515 grad_grads[start + n_curls][0][2][0] =
516 -monoval_jplus[0][1] * monoval_i[1][1] * monoval[2][1];
517 grad_grads[start + n_curls][0][2][1] =
518 -monoval_jplus[0][0] * monoval_i[1][2] * monoval[2][1];
519 grad_grads[start + n_curls][0][2][2] =
520 -monoval_jplus[0][0] * monoval_i[1][1] * monoval[2][2];
521 grad_grads[start + n_curls][1][0][0] =
522 monoval_j[0][2] * monoval_i[1][0] * monoval[2][0] *
523 static_cast<double>(j + my_degree + 2);
524 grad_grads[start + n_curls][1][0][1] =
525 monoval_j[0][1] * monoval_i[1][1] * monoval[2][0] *
526 static_cast<double>(j + my_degree + 2);
527 grad_grads[start + n_curls][1][0][2] =
528 monoval_j[0][1] * monoval_i[1][0] * monoval[2][1] *
529 static_cast<double>(j + my_degree + 2);
530 grad_grads[start + n_curls][1][1][0] =
531 monoval_j[0][1] * monoval_i[1][1] * monoval[2][0] *
532 static_cast<double>(j + my_degree + 2);
533 grad_grads[start + n_curls][1][1][1] =
534 monoval_j[0][0] * monoval_i[1][2] * monoval[2][0] *
535 static_cast<double>(j + my_degree + 2);
536 grad_grads[start + n_curls][1][1][2] =
537 monoval_j[0][0] * monoval_i[1][1] * monoval[2][1] *
538 static_cast<double>(j + my_degree + 2);
539 grad_grads[start + n_curls][1][2][0] =
540 monoval_j[0][1] * monoval_i[1][0] * monoval[2][1] *
541 static_cast<double>(j + my_degree + 2);
542 grad_grads[start + n_curls][1][2][1] =
543 monoval_j[0][0] * monoval_i[1][1] * monoval[2][1] *
544 static_cast<double>(j + my_degree + 2);
545 grad_grads[start + n_curls][1][2][2] =
546 monoval_j[0][0] * monoval_i[1][0] * monoval[2][2] *
547 static_cast<double>(j + my_degree + 2);
548 grad_grads[start + n_curls][2][0][0] =
549 -monoval_j[0][2] * monoval_i[1][1] * monoval_plus[2][0];
550 grad_grads[start + n_curls][2][0][1] =
551 -monoval_j[0][1] * monoval_i[1][2] * monoval_plus[2][0];
552 grad_grads[start + n_curls][2][0][2] =
553 -monoval_j[0][1] * monoval_i[1][1] * monoval_plus[2][1];
554 grad_grads[start + n_curls][2][1][0] =
555 -monoval_j[0][1] * monoval_i[1][2] * monoval_plus[2][0];
556 grad_grads[start + n_curls][2][1][1] =
557 -monoval_j[0][0] * monoval_i[1][3] * monoval_plus[2][0];
558 grad_grads[start + n_curls][2][1][2] =
559 -monoval_j[0][0] * monoval_i[1][2] * monoval_plus[2][1];
560 grad_grads[start + n_curls][2][2][0] =
561 -monoval_j[0][1] * monoval_i[1][1] * monoval_plus[2][1];
562 grad_grads[start + n_curls][2][2][1] =
563 -monoval_j[0][0] * monoval_i[1][2] * monoval_plus[2][1];
564 grad_grads[start + n_curls][2][2][2] =
565 -monoval_j[0][0] * monoval_i[1][1] * monoval_plus[2][2];
567 grad_grads[start + 2 * n_curls][0][0][0] =
568 -monoval_jplus[0][2] * monoval[1][0] * monoval_i[2][1];
569 grad_grads[start + 2 * n_curls][0][0][1] =
570 -monoval_jplus[0][1] * monoval[1][1] * monoval_i[2][1];
571 grad_grads[start + 2 * n_curls][0][0][2] =
572 -monoval_jplus[0][1] * monoval[1][0] * monoval_i[2][2];
573 grad_grads[start + 2 * n_curls][0][1][0] =
574 -monoval_jplus[0][1] * monoval[1][1] * monoval_i[2][1];
575 grad_grads[start + 2 * n_curls][0][1][1] =
576 -monoval_jplus[0][0] * monoval[1][2] * monoval_i[2][1];
577 grad_grads[start + 2 * n_curls][0][1][2] =
578 -monoval_jplus[0][0] * monoval[1][1] * monoval_i[2][2];
579 grad_grads[start + 2 * n_curls][0][2][0] =
580 -monoval_jplus[0][1] * monoval[1][0] * monoval_i[2][2];
581 grad_grads[start + 2 * n_curls][0][2][1] =
582 -monoval_jplus[0][0] * monoval[1][1] * monoval_i[2][2];
583 grad_grads[start + 2 * n_curls][0][2][2] =
584 -monoval_jplus[0][0] * monoval[1][0] * monoval_i[2][3];
585 grad_grads[start + 2 * n_curls][1][0][0] =
586 -monoval_j[0][2] * monoval_plus[1][0] * monoval_i[2][1];
587 grad_grads[start + 2 * n_curls][1][0][1] =
588 -monoval_j[0][1] * monoval_plus[1][1] * monoval_i[2][1];
589 grad_grads[start + 2 * n_curls][1][0][2] =
590 -monoval_j[0][1] * monoval_plus[1][0] * monoval_i[2][2];
591 grad_grads[start + 2 * n_curls][1][1][0] =
592 -monoval_j[0][1] * monoval_plus[1][1] * monoval_i[2][1];
593 grad_grads[start + 2 * n_curls][1][1][1] =
594 -monoval_j[0][0] * monoval_plus[1][2] * monoval_i[2][1];
595 grad_grads[start + 2 * n_curls][1][1][2] =
596 -monoval_j[0][0] * monoval_plus[1][1] * monoval_i[2][2];
597 grad_grads[start + 2 * n_curls][1][2][0] =
598 -monoval_j[0][1] * monoval_plus[1][0] * monoval_i[2][2];
599 grad_grads[start + 2 * n_curls][1][2][1] =
600 -monoval_j[0][0] * monoval_plus[1][1] * monoval_i[2][2];
601 grad_grads[start + 2 * n_curls][1][2][2] =
602 -monoval_j[0][0] * monoval_plus[1][0] * monoval_i[2][3];
603 grad_grads[start + 2 * n_curls][2][0][0] =
604 monoval_j[0][2] * monoval[1][0] * monoval_i[2][0] *
605 static_cast<double>(j + my_degree + 2);
606 grad_grads[start + 2 * n_curls][2][0][1] =
607 monoval_j[0][1] * monoval[1][1] * monoval_i[2][0] *
608 static_cast<double>(j + my_degree + 2);
609 grad_grads[start + 2 * n_curls][2][0][2] =
610 monoval_j[0][1] * monoval[1][0] * monoval_i[2][1] *
611 static_cast<double>(j + my_degree + 2);
612 grad_grads[start + 2 * n_curls][2][1][0] =
613 monoval_j[0][1] * monoval[1][1] * monoval_i[2][0] *
614 static_cast<double>(j + my_degree + 2);
615 grad_grads[start + 2 * n_curls][2][1][1] =
616 monoval_j[0][0] * monoval[1][2] * monoval_i[2][0] *
617 static_cast<double>(j + my_degree + 2);
618 grad_grads[start + 2 * n_curls][2][1][2] =
619 monoval_j[0][0] * monoval[1][1] * monoval_i[2][1] *
620 static_cast<double>(j + my_degree + 2);
621 grad_grads[start + 2 * n_curls][2][2][0] =
622 monoval_j[0][1] * monoval[1][0] * monoval_i[2][1] *
623 static_cast<double>(j + my_degree + 2);
624 grad_grads[start + 2 * n_curls][2][2][1] =
625 monoval_j[0][0] * monoval[1][1] * monoval_i[2][1] *
626 static_cast<double>(j + my_degree + 2);
627 grad_grads[start + 2 * n_curls][2][2][2] =
628 monoval_j[0][0] * monoval[1][0] * monoval_i[2][2] *
629 static_cast<double>(j + my_degree + 2);
633 grad_grads[start + 1][0][0][0] =
634 monoval_i[0][2] * monoval[1][0] * monoval_j[2][0] *
635 static_cast<double>(j + my_degree + 2);
636 grad_grads[start + 1][0][0][1] =
637 monoval_i[0][1] * monoval[1][1] * monoval_j[2][0] *
638 static_cast<double>(j + my_degree + 2);
639 grad_grads[start + 1][0][0][2] =
640 monoval_i[0][1] * monoval[1][0] * monoval_j[2][1] *
641 static_cast<double>(j + my_degree + 2);
642 grad_grads[start + 1][0][1][0] =
643 monoval_i[0][1] * monoval[1][1] * monoval_j[2][0] *
644 static_cast<double>(j + my_degree + 2);
645 grad_grads[start + 1][0][1][1] =
646 monoval_i[0][0] * monoval[1][2] * monoval_j[2][0] *
647 static_cast<double>(j + my_degree + 2);
648 grad_grads[start + 1][0][1][2] =
649 monoval_i[0][0] * monoval[1][1] * monoval_j[2][1] *
650 static_cast<double>(j + my_degree + 2);
651 grad_grads[start + 1][0][2][0] =
652 monoval_i[0][1] * monoval[1][0] * monoval_j[2][1] *
653 static_cast<double>(j + my_degree + 2);
654 grad_grads[start + 1][0][2][1] =
655 monoval_i[0][0] * monoval[1][1] * monoval_j[2][1] *
656 static_cast<double>(j + my_degree + 2);
657 grad_grads[start + 1][0][2][2] =
658 monoval_i[0][0] * monoval[1][0] * monoval_j[2][2] *
659 static_cast<double>(j + my_degree + 2);
660 grad_grads[start + 1][1][0][0] =
661 -monoval_i[0][3] * monoval_plus[1][0] * monoval_j[2][0];
662 grad_grads[start + 1][1][0][1] =
663 -monoval_i[0][2] * monoval_plus[1][1] * monoval_j[2][0];
664 grad_grads[start + 1][1][0][2] =
665 -monoval_i[0][2] * monoval_plus[1][0] * monoval_j[2][1];
666 grad_grads[start + 1][1][1][0] =
667 -monoval_i[0][2] * monoval_plus[1][1] * monoval_j[2][0];
668 grad_grads[start + 1][1][1][1] =
669 -monoval_i[0][1] * monoval_plus[1][2] * monoval_j[2][0];
670 grad_grads[start + 1][1][1][2] =
671 -monoval_i[0][1] * monoval_plus[1][1] * monoval_j[2][1];
672 grad_grads[start + 1][1][2][0] =
673 -monoval_i[0][2] * monoval_plus[1][0] * monoval_j[2][1];
674 grad_grads[start + 1][1][2][1] =
675 -monoval_i[0][1] * monoval_plus[1][1] * monoval_j[2][1];
676 grad_grads[start + 1][1][2][2] =
677 -monoval_i[0][1] * monoval_plus[1][0] * monoval_j[2][2];
678 grad_grads[start + 1][2][0][0] =
679 -monoval_i[0][3] * monoval[1][0] * monoval_jplus[2][0];
680 grad_grads[start + 1][2][0][1] =
681 -monoval_i[0][2] * monoval[1][1] * monoval_jplus[2][0];
682 grad_grads[start + 1][2][0][2] =
683 -monoval_i[0][2] * monoval[1][0] * monoval_jplus[2][1];
684 grad_grads[start + 1][2][1][0] =
685 -monoval_i[0][2] * monoval[1][1] * monoval_jplus[2][0];
686 grad_grads[start + 1][2][1][1] =
687 -monoval_i[0][1] * monoval[1][2] * monoval_jplus[2][0];
688 grad_grads[start + 1][2][1][2] =
689 -monoval_i[0][1] * monoval[1][1] * monoval_jplus[2][1];
690 grad_grads[start + 1][2][2][0] =
691 -monoval_i[0][2] * monoval[1][0] * monoval_jplus[2][1];
692 grad_grads[start + 1][2][2][1] =
693 -monoval_i[0][1] * monoval[1][1] * monoval_jplus[2][1];
694 grad_grads[start + 1][2][2][2] =
695 -monoval_i[0][1] * monoval[1][0] * monoval_jplus[2][2];
697 grad_grads[start + n_curls + 1][0][0][0] =
698 -monoval_plus[0][2] * monoval_i[1][1] * monoval_j[2][0];
699 grad_grads[start + n_curls + 1][0][0][1] =
700 -monoval_plus[0][1] * monoval_i[1][2] * monoval_j[2][0];
701 grad_grads[start + n_curls + 1][0][0][2] =
702 -monoval_plus[0][1] * monoval_i[1][1] * monoval_j[2][1];
703 grad_grads[start + n_curls + 1][0][1][0] =
704 -monoval_plus[0][1] * monoval_i[1][2] * monoval_j[2][0];
705 grad_grads[start + n_curls + 1][0][1][1] =
706 -monoval_plus[0][0] * monoval_i[1][3] * monoval_j[2][0];
707 grad_grads[start + n_curls + 1][0][1][2] =
708 -monoval_plus[0][0] * monoval_i[1][2] * monoval_j[2][1];
709 grad_grads[start + n_curls + 1][0][2][0] =
710 -monoval_plus[0][1] * monoval_i[1][1] * monoval_j[2][1];
711 grad_grads[start + n_curls + 1][0][2][1] =
712 -monoval_plus[0][0] * monoval_i[1][2] * monoval_j[2][1];
713 grad_grads[start + n_curls + 1][0][2][2] =
714 -monoval_plus[0][0] * monoval_i[1][1] * monoval_j[2][2];
715 grad_grads[start + n_curls + 1][1][0][0] =
716 monoval[0][2] * monoval_i[1][0] * monoval_j[2][0] *
717 static_cast<double>(j + my_degree + 2);
718 grad_grads[start + n_curls + 1][1][0][1] =
719 monoval[0][1] * monoval_i[1][1] * monoval_j[2][0] *
720 static_cast<double>(j + my_degree + 2);
721 grad_grads[start + n_curls + 1][1][0][2] =
722 monoval[0][1] * monoval_i[1][0] * monoval_j[2][1] *
723 static_cast<double>(j + my_degree + 2);
724 grad_grads[start + n_curls + 1][1][1][0] =
725 monoval[0][1] * monoval_i[1][1] * monoval_j[2][0] *
726 static_cast<double>(j + my_degree + 2);
727 grad_grads[start + n_curls + 1][1][1][1] =
728 monoval[0][0] * monoval_i[1][2] * monoval_j[2][0] *
729 static_cast<double>(j + my_degree + 2);
730 grad_grads[start + n_curls + 1][1][1][2] =
731 monoval[0][0] * monoval_i[1][1] * monoval_j[2][1] *
732 static_cast<double>(j + my_degree + 2);
733 grad_grads[start + n_curls + 1][1][2][0] =
734 monoval[0][1] * monoval_i[1][0] * monoval_j[2][1] *
735 static_cast<double>(j + my_degree + 2);
736 grad_grads[start + n_curls + 1][1][2][1] =
737 monoval[0][0] * monoval_i[1][1] * monoval_j[2][1] *
738 static_cast<double>(j + my_degree + 2);
739 grad_grads[start + n_curls + 1][1][2][2] =
740 monoval[0][0] * monoval_i[1][0] * monoval_j[2][2] *
741 static_cast<double>(j + my_degree + 2);
742 grad_grads[start + n_curls + 1][2][0][0] =
743 -monoval[0][2] * monoval_i[1][1] * monoval_jplus[2][0];
744 grad_grads[start + n_curls + 1][2][0][1] =
745 -monoval[0][1] * monoval_i[1][2] * monoval_jplus[2][0];
746 grad_grads[start + n_curls + 1][2][0][2] =
747 -monoval[0][1] * monoval_i[1][1] * monoval_jplus[2][1];
748 grad_grads[start + n_curls + 1][2][1][0] =
749 -monoval[0][1] * monoval_i[1][2] * monoval_jplus[2][0];
750 grad_grads[start + n_curls + 1][2][1][1] =
751 -monoval[0][0] * monoval_i[1][3] * monoval_jplus[2][0];
752 grad_grads[start + n_curls + 1][2][1][2] =
753 -monoval[0][0] * monoval_i[1][2] * monoval_jplus[2][1];
754 grad_grads[start + n_curls + 1][2][2][0] =
755 -monoval[0][1] * monoval_i[1][1] * monoval_jplus[2][1];
756 grad_grads[start + n_curls + 1][2][2][1] =
757 -monoval[0][0] * monoval_i[1][2] * monoval_jplus[2][1];
758 grad_grads[start + n_curls + 1][2][2][2] =
759 -monoval[0][0] * monoval_i[1][1] * monoval_jplus[2][2];
761 grad_grads[start + 2 * n_curls + 1][0][0][0] =
762 -monoval_plus[0][2] * monoval_j[1][0] * monoval_i[2][1];
763 grad_grads[start + 2 * n_curls + 1][0][0][1] =
764 -monoval_plus[0][1] * monoval_j[1][1] * monoval_i[2][1];
765 grad_grads[start + 2 * n_curls + 1][0][0][2] =
766 -monoval_plus[0][1] * monoval_j[1][0] * monoval_i[2][2];
767 grad_grads[start + 2 * n_curls + 1][0][1][0] =
768 -monoval_plus[0][1] * monoval_j[1][1] * monoval_i[2][1];
769 grad_grads[start + 2 * n_curls + 1][0][1][1] =
770 -monoval_plus[0][0] * monoval_j[1][2] * monoval_i[2][1];
771 grad_grads[start + 2 * n_curls + 1][0][1][2] =
772 -monoval_plus[0][0] * monoval_j[1][1] * monoval_i[2][2];
773 grad_grads[start + 2 * n_curls + 1][0][2][0] =
774 -monoval_plus[0][1] * monoval_j[1][0] * monoval_i[2][2];
775 grad_grads[start + 2 * n_curls + 1][0][2][1] =
776 -monoval_plus[0][0] * monoval_j[1][1] * monoval_i[2][2];
777 grad_grads[start + 2 * n_curls + 1][0][2][2] =
778 -monoval_plus[0][0] * monoval_j[1][0] * monoval_i[2][3];
779 grad_grads[start + 2 * n_curls + 1][1][0][0] =
780 -monoval[0][2] * monoval_jplus[1][0] * monoval_i[2][1];
781 grad_grads[start + 2 * n_curls + 1][1][0][1] =
782 -monoval[0][1] * monoval_jplus[1][1] * monoval_i[2][1];
783 grad_grads[start + 2 * n_curls + 1][1][0][2] =
784 -monoval[0][1] * monoval_jplus[1][0] * monoval_i[2][2];
785 grad_grads[start + 2 * n_curls + 1][1][1][0] =
786 -monoval[0][1] * monoval_jplus[1][1] * monoval_i[2][1];
787 grad_grads[start + 2 * n_curls + 1][1][1][1] =
788 -monoval[0][0] * monoval_jplus[1][2] * monoval_i[2][1];
789 grad_grads[start + 2 * n_curls + 1][1][1][2] =
790 -monoval[0][0] * monoval_jplus[1][1] * monoval_i[2][2];
791 grad_grads[start + 2 * n_curls + 1][1][2][0] =
792 -monoval[0][1] * monoval_jplus[1][0] * monoval_i[2][2];
793 grad_grads[start + 2 * n_curls + 1][1][2][1] =
794 -monoval[0][0] * monoval_jplus[1][1] * monoval_i[2][2];
795 grad_grads[start + 2 * n_curls + 1][1][2][2] =
796 -monoval[0][0] * monoval_jplus[1][0] * monoval_i[2][3];
797 grad_grads[start + 2 * n_curls + 1][2][0][0] =
798 monoval[0][2] * monoval_j[1][0] * monoval_i[2][0] *
799 static_cast<double>(j + my_degree + 2);
800 grad_grads[start + 2 * n_curls + 1][2][0][1] =
801 monoval[0][1] * monoval_j[1][1] * monoval_i[2][0] *
802 static_cast<double>(j + my_degree + 2);
803 grad_grads[start + 2 * n_curls + 1][2][0][2] =
804 monoval[0][1] * monoval_j[1][0] * monoval_i[2][1] *
805 static_cast<double>(j + my_degree + 2);
806 grad_grads[start + 2 * n_curls + 1][2][1][0] =
807 monoval[0][1] * monoval_j[1][1] * monoval_i[2][0] *
808 static_cast<double>(j + my_degree + 2);
809 grad_grads[start + 2 * n_curls + 1][2][1][1] =
810 monoval[0][0] * monoval_j[1][2] * monoval_i[2][0] *
811 static_cast<double>(j + my_degree + 2);
812 grad_grads[start + 2 * n_curls + 1][2][1][2] =
813 monoval[0][0] * monoval_j[1][1] * monoval_i[2][1] *
814 static_cast<double>(j + my_degree + 2);
815 grad_grads[start + 2 * n_curls + 1][2][2][0] =
816 monoval[0][1] * monoval_j[1][0] * monoval_i[2][1] *
817 static_cast<double>(j + my_degree + 2);
818 grad_grads[start + 2 * n_curls + 1][2][2][1] =
819 monoval[0][0] * monoval_j[1][1] * monoval_i[2][1] *
820 static_cast<double>(j + my_degree + 2);
821 grad_grads[start + 2 * n_curls + 1][2][2][2] =
822 monoval[0][0] * monoval_j[1][0] * monoval_i[2][2] *
823 static_cast<double>(j + my_degree + 2);