proteus 1.9.0
C/C++/Fortran libraries
Loading...
Searching...
No Matches
ADR.h
Go to the documentation of this file.
1#ifndef ADR_H
2#define ADR_H
3#include <cmath>
4#include <iostream>
5#include <set>
6#include <stdexcept>
7#include <map>
8#include <valarray>
9#include <vector>
10#include "CompKernel.h"
11#include "ModelFactory.h"
12#include "equivalent_polynomials.h"
13#include "mprans/ArgumentsDict.h"
14#include "xtensor-python/pyarray.hpp"
15
16namespace py = pybind11;
17
18namespace proteus
19{
20 template <int nSpace, int nP_ifem, int nP, int nQ, int nEBQ>
21 using GeneralizedFunctions = equivalent_polynomials::GeneralizedFunctions_mix<nSpace, nP_ifem, nP, nQ, nEBQ, true>;
22
23 class cADR_base
24 {
25 public:
26 virtual ~cADR_base() {}
27 virtual void calculateResidual(arguments_dict &args) = 0;
28 virtual void calculateJacobian(arguments_dict &args) = 0;
29 };
30
31 template <class CompKernelType,
32 int nSpace,
33 int nQuadraturePoints_element,
34 int nDOF_mesh_trial_element,
35 int nDOF_trial_element,
36 int nDOF_test_element,
37 int nQuadraturePoints_elementBoundary>
38 class cADR : public cADR_base
39 {
40 public:
41 std::set<int> ifem_boundaries, ifem_boundary_elements,
42 cutfem_boundaries, cutfem_boundary_elements;
43 std::valarray<bool> elementIsActive;
44 const int nDOF_test_X_trial_element;
45 CompKernelType ck;
46 // nEBQ is sized for ALL faces (nDOF_mesh_trial_element * nQuadraturePoints_elementBoundary),
47 // not one face, so the immersed-interface facet terms can evaluate va/vb at the quadrature
48 // points of the face they are actually integrating -- see the xB_ref_faces comment in
49 // calculateResidual for why passing the reference *boundary* points directly is wrong.
50 using GfType = GeneralizedFunctions<nSpace, nDOF_trial_element, 4, nQuadraturePoints_element, nDOF_mesh_trial_element * nQuadraturePoints_elementBoundary>;
51 // Degree of the edge-restricted moment-fit Heaviside used by the immersed-interface facet
52 // terms; must match GfType's nP above so the fit covers the facet integrands exactly
53 // (degree <= 1 for P1, <= 3 for P2).
54 static const int nP_edge = 4;
55 // Per-element cache of the equivalent-polynomial/IFEM reconstruction
56 // (gf_s.calculate()/gf_f.calculate()): permutation, cut classification,
57 // the IFEM basis coefficient solve, and the H/ImH/D + VA/VB (+
58 // gradients) evaluated at every quadrature point. None of this depends
59 // on the current solution u -- only on the level-set geometry
60 // (element_phi_s/element_phi_f) and mua/mub/jf -- so it is only
61 // recomputed when ifemGeometryGeneration advances (see
62 // recomputeIFEMGeometry / markIFEMGeometryDirty on the Python side).
63 // Interior (isBoundary=false) and boundary (isBoundary=true) evaluations
64 // populate different internal storage within the same GfType object, so
65 // they are tracked with separate generation/icase arrays.
66 std::vector<GfType> gf_f_cache, gf_s_cache;
67 std::vector<int> gf_s_interior_gen, gf_f_interior_gen, gf_f_boundary_gen;
68 std::vector<int> gf_s_interior_icase, gf_f_interior_icase, gf_f_boundary_icase;
69 int ifemGeometryGeneration = 0;
70 cADR() : nDOF_test_X_trial_element(nDOF_test_element * nDOF_trial_element), ck()
71 {
72 }
73 inline void ensureIFEMCacheSized(int nElements_global)
74 {
75 if ((int)gf_f_cache.size() == nElements_global)
76 return;
77 gf_f_cache.assign(nElements_global, GfType());
78 gf_s_cache.assign(nElements_global, GfType());
79 gf_s_interior_gen.assign(nElements_global, -1);
80 gf_f_interior_gen.assign(nElements_global, -1);
81 gf_f_boundary_gen.assign(nElements_global, -1);
82 gf_s_interior_icase.assign(nElements_global, 0);
83 gf_f_interior_icase.assign(nElements_global, 0);
84 gf_f_boundary_icase.assign(nElements_global, 0);
85 for (auto &gf : gf_f_cache)
86 gf.useExact = true;
87 for (auto &gf : gf_s_cache)
88 gf.useExact = true;
89 }
90
91 inline void exteriorNumericalDiffusiveFlux(int *rowptr,
92 int *colind,
93 const int &isDOFBoundary,
94 const int &isDiffusiveFluxBoundary,
95 const double n[nSpace],
96 double *bc_a,
97 const double &bc_u,
98 const double &bc_flux,
99 double *a,
100 const double grad_potential[nSpace],
101 const double &u,
102 const double &penalty,
103 double &flux)
104 {
105 double diffusiveVelocityComponent_I;
106 double penaltyFlux;
107 double max_a;
108 if (isDiffusiveFluxBoundary == 1)
109 {
110 flux = bc_flux;
111 }
112 else if (isDOFBoundary == 1)
113 {
114 flux = 0.0;
115 max_a = 0.0;
116 for (int I = 0; I < nSpace; I++)
117 {
118 diffusiveVelocityComponent_I = 0.0;
119 for (int m = rowptr[I]; m < rowptr[I + 1]; m++)
120 {
121 diffusiveVelocityComponent_I -= a[m] * grad_potential[colind[m]];
122 max_a = fmax(max_a, a[m]);
123 }
124 flux += diffusiveVelocityComponent_I * n[I];
125 }
126 penaltyFlux = max_a * penalty * (u - bc_u);
127 flux += penaltyFlux;
128 }
129 else
130 {
131 std::cerr << "warning, diffusion term with no boundary condition set, setting diffusive flux to 0.0" << std::endl;
132 flux = 0.0;
133 }
134 }
135
136 inline double ExteriorNumericalDiffusiveFluxJacobian(int *rowptr,
137 int *colind,
138 const int &isDOFBoundary,
139 const int &isDiffusiveFluxBoundary,
140 const double n[nSpace],
141 double *a,
142 const double &v,
143 const double grad_v[nSpace],
144 const double &penalty)
145 {
146 double dvel_I;
147 double tmp = 0.0;
148 double max_a = 0.0;
149 if ((isDiffusiveFluxBoundary == 0) && (isDOFBoundary == 1))
150 {
151 for (int I = 0; I < nSpace; I++)
152 {
153 dvel_I = 0.0;
154 for (int m = rowptr[I]; m < rowptr[I + 1]; m++)
155 {
156 dvel_I -= a[m] * grad_v[colind[m]];
157 max_a = fmax(max_a, a[m]);
158 }
159 tmp += dvel_I * n[I];
160 }
161 tmp += max_a * penalty * v;
162 }
163 return tmp;
164 }
165
166 inline void calculateSubgridError_tau(const double &elementDiameter,
167 const double &dmt,
168 const double dH[nSpace],
169 double &cfl,
170 double &tau)
171 {
172 double h;
173 double nrm_v;
174 double oneByAbsdt;
175 h = elementDiameter;
176 nrm_v = 0.0;
177 for (int I = 0; I < nSpace; I++)
178 nrm_v += dH[I] * dH[I];
179 nrm_v = sqrt(nrm_v);
180 cfl = nrm_v / h;
181 oneByAbsdt = fabs(dmt);
182 tau = 1.0 / (2.0 * nrm_v / h + oneByAbsdt + 1.0e-8);
183 }
184
185 inline void calculateSubgridError_tau(const double &Ct_sge,
186 const double G[nSpace * nSpace],
187 const double &A0,
188 const double Ai[nSpace],
189 double &tau_v,
190 double &cfl)
191 {
192 double v_d_Gv = 0.0;
193 for (int I = 0; I < nSpace; I++)
194 for (int J = 0; J < nSpace; J++)
195 v_d_Gv += Ai[I] * G[I * nSpace + J] * Ai[J];
196 tau_v = 1.0 / sqrt(Ct_sge * A0 * A0 + v_d_Gv + 1.0e-8);
197 }
198
199 inline void calculateNumericalDiffusion(const double &shockCapturingDiffusion,
200 const double &elementDiameter,
201 const double &strong_residual,
202 const double grad_u[nSpace],
203 double &numDiff)
204 {
205 double h;
206 double num;
207 double den;
208 double n_grad_u;
209 h = elementDiameter;
210 n_grad_u = 0.0;
211 for (int I = 0; I < nSpace; I++)
212 n_grad_u += grad_u[I] * grad_u[I];
213 num = shockCapturingDiffusion * 0.5 * h * fabs(strong_residual);
214 den = sqrt(n_grad_u) + 1.0e-8;
215 numDiff = num / den;
216 }
217
218 inline void exteriorNumericalAdvectiveFlux(const int &isDOFBoundary_u,
219 const int &isFluxBoundary_u,
220 const double n[nSpace],
221 const double &bc_u,
222 const double &bc_flux_u,
223 const double &u,
224 const double velocity[nSpace],
225 double &flux)
226 {
227
228 double flow = 0.0;
229 for (int I = 0; I < nSpace; I++)
230 flow += n[I] * velocity[I];
231 if (isDOFBoundary_u == 1)
232 {
233 if (flow >= 0.0)
234 {
235 flux = u * flow;
236 // flux = flow;
237 }
238 else
239 {
240 flux = bc_u * flow;
241 // flux = flow;
242 }
243 }
244 else if (isFluxBoundary_u == 1)
245 {
246 flux = bc_flux_u;
247 }
248 else
249 {
250 if (flow >= 0.0)
251 {
252 flux = u * flow;
253 }
254 else
255 {
256 flux = 0.0;
257 }
258 }
259 // flux = flow;
260 }
261
262 inline void exteriorNumericalAdvectiveFluxDerivative(const int &isDOFBoundary_u,
263 const int &isFluxBoundary_u,
264 const double n[nSpace],
265 const double velocity[nSpace],
266 double &dflux)
267 {
268 double flow = 0.0;
269 for (int I = 0; I < nSpace; I++)
270 {
271 flow += n[I] * velocity[I];
272 }
273 // double flow=n[0]*velocity[0]+n[1]*velocity[1]+n[2]*velocity[2];
274 dflux = 0.0; // default to no flux
275 if (isDOFBoundary_u == 1)
276 {
277 if (flow >= 0.0)
278 {
279 dflux = flow;
280 }
281 else
282 {
283 dflux = 0.0;
284 }
285 }
286 else if (isFluxBoundary_u == 1)
287 {
288 dflux = 0.0;
289 }
290 else
291 {
292 if (flow >= 0.0)
293 {
294 dflux = flow;
295 }
296 }
297 }
298
299 inline void updateEmbeddedBoundaryTerms(const double embeddedBoundary_penalty,
300 const double dV,
301 double *embeddedBoundary_normal,
302 const double u_s,
303 const double u,
304 const double grad_u[nSpace],
305 const double a,
306 double &r,
307 double &dr,
308 double &ham,
309 double *dham,
310 double *f,
311 double *df,
312 const double D_s)
313 {
314 // todo this doesn't look 1d/3d
315 double outward_normal[nSpace];
316 for (int I = 0; I < nSpace; I++)
317 outward_normal[I] = -embeddedBoundary_normal[I];
318
319 // diffusive flux
320 for (int I = 0; I < nSpace; I++)
321 {
322 ham -= outward_normal[I] * grad_u[I];
323 dham[I] -= D_s * a * outward_normal[I];
324 // Nitsche adjoint consistency
325 f[I] += D_s * a * outward_normal[I] * (u - u_s);
326 df[I] += D_s * a * outward_normal[I];
327 }
328 ham *= D_s * a;
329
330 // Nitsche Dirichlet penalty
331 r += D_s * a * embeddedBoundary_penalty * (u - u_s);
332 dr += D_s * a * embeddedBoundary_penalty;
333 }
334
335 inline void updateImmersedBoundaryTerms(const double immersedBoundary_penalty,
336 const double dV,
337 double *immersedBoundary_normal,
338 double x,
339 double y,
340 double z,
341 const double u_s,
342 const double u,
343 const double grad_u[nSpace],
344 const double a,
345 double &r,
346 double &dr,
347 double &ham,
348 double *dham,
349 double *f,
350 double *df,
351 const double fluxJump,
352 const double *fluxJumpVector,
353 const double D_f)
354 {
355 // todo this doesn't look 1d/3d
356 /* double outward_normal[nSpace];
357 for (int I=0;I<nSpace;I++)
358 outward_normal[I] = immersedBoundary_normal[I];
359
360
361
362 //diffusive flux
363 for (int I=0;I<nSpace;I++)
364 {
365 ham -= outward_normal[I] * grad_u[I];
366 dham[I] -= D_f * a * outward_normal[I];
367 //Nitsche adjoint consistency
368 f[I] += D_f * a * outward_normal[I] * (u - u_s);
369 df[I] += D_f * a * outward_normal[I];
370 }
371 ham *= D_f * a;
372 //Nitsche Dirichlet penalty
373 r += D_f*a*immersedBoundary_penalty * (u - u_s);
374 dr += D_f*a*immersedBoundary_penalty; */
375 // std::cout << "D_f = " << D_f << std::endl;
376 // Prescribed flux jump across the interface, supplied by the physics file as
377 // [beta du/dn] = fluxJump(x) + fluxJumpVector(x) . n. The vector part carries
378 // the cases whose jump depends on the interface normal. Both default to zero,
379 // so a problem with no prescribed flux jump contributes nothing here.
380 double jump_flux = fluxJump;
381 for (int I = 0; I < nSpace; I++)
382 jump_flux += fluxJumpVector[I] * immersedBoundary_normal[I];
383 r += jump_flux * D_f;
384 dr = 0.0;
385 ham = 0.0;
386 dham[0] = 0.0;
387 dham[1] = 0.0;
388 f[0] = 0.0;
389 f[1] = 0.0;
390 df[0] = 0.0;
391 df[1] = 0.0;
392 }
393
394 inline void calculateElementResidual(int icase_f,
395 // element
396 xt::pyarray<double> &mesh_trial_ref,
397 xt::pyarray<double> &mesh_grad_trial_ref,
398 xt::pyarray<double> &mesh_dof,
399 xt::pyarray<int> &mesh_l2g,
400 xt::pyarray<double> &x_ref,
401 xt::pyarray<double> &dV_ref,
402 xt::pyarray<double> &u_trial_ref,
403 xt::pyarray<double> &u_grad_trial_ref,
404 xt::pyarray<double> &u_test_ref,
405 xt::pyarray<double> &u_grad_test_ref,
406 xt::pyarray<double> &elementDiameter,
407 xt::pyarray<double> &elementBoundaryDiameter,
408 xt::pyarray<double> &nodeDiametersArray,
409 xt::pyarray<double> &cfl,
410 double Ct_sge,
411 double sc_uref,
412 double sc_alpha,
413 double useMetrics,
414 // element boundary
415 xt::pyarray<double> &mesh_trial_trace_ref,
416 xt::pyarray<double> &mesh_grad_trial_trace_ref,
417 xt::pyarray<double> &dS_ref,
418 xt::pyarray<double> &u_trial_trace_ref,
419 xt::pyarray<double> &u_grad_trial_trace_ref,
420 xt::pyarray<double> &u_test_trace_ref,
421 xt::pyarray<double> &u_grad_test_trace_ref,
422 xt::pyarray<double> &normal_ref,
423 xt::pyarray<double> &boundaryJac_ref,
424 // physics
425 int nElements_global,
426 int nElementBoundaries_owned,
427 xt::pyarray<int> &u_l2g,
428 xt::pyarray<double> &u_dof,
429 xt::pyarray<int> &sd_rowptr,
430 xt::pyarray<int> &sd_colind,
431 xt::pyarray<double> &q_a,
432 xt::pyarray<double> &q_v,
433 xt::pyarray<double> &q_r,
434 int lag_shockCapturingDiffusion,
435 double shockCapturingDiffusion,
436 xt::pyarray<double> &q_numDiff_u,
437 xt::pyarray<double> &q_numDiff_u_last,
438 int offset_u,
439 int stride_u,
440 xt::pyarray<double> &elementResidual_u,
441 int nExteriorElementBoundaries_global,
442 xt::pyarray<int> &exteriorElementBoundariesArray,
443 xt::pyarray<int> &elementBoundariesArray,
444 xt::pyarray<int> &elementBoundaryElementsArray,
445 xt::pyarray<int> &elementBoundaryLocalElementBoundariesArray,
446 xt::pyarray<double> &element_u,
447 int eN,
448 const bool embeddedBoundary,
449 const double embeddedBoundary_penalty,
450 xt::pyarray<double> &embeddedBoundary_normal_q,
451 xt::pyarray<double> &embeddedBoundary_u_q,
452 const bool immersedBoundary,
453 const double immersedBoundary_penalty,
454 xt::pyarray<double> &immersedBoundary_sdf_q,
455 xt::pyarray<double> &immersedBoundary_normal_q,
456 xt::pyarray<double> &immersedBoundary_u_q,
457 xt::pyarray<double> &immersedBoundary_fluxJump_q,
458 xt::pyarray<double> &immersedBoundary_fluxJumpVector_q,
459 xt::pyarray<double> &immersedBoundary_solutionJump_nodes,
460 double *element_phi_f,
461 bool &element_active,
462 std::valarray<bool> &elementIsActive,
463 double *JA,
464 double *JB,
465 double &L2_error,
466 double &Linfty_error,
467 double mua,
468 double mub,
469 xt::pyarray<double> &q_u_exact_inner,
470 xt::pyarray<double> &q_u_exact_outer,
471 const bool PG)
472 {
473 // per-element cached equivalent-polynomial/IFEM reconstruction
474 // (see ensureIFEMCacheSized / ifemGeometryGeneration)
475 GfType &gf_f = gf_f_cache[eN];
476 GfType &gf_s = gf_s_cache[eN];
477 for (int i = 0; i < nDOF_test_element; i++)
478 {
479 elementResidual_u.data()[i] = 0.0;
480 }
481 // std::cout << "Calculating element residual for element " << eN << std::endl;
482 // loop over quadrature points and compute integrands
483 for (int k = 0; k < nQuadraturePoints_element; k++)
484 {
485 // std::cout << " quadrature point " << k << "\t" << x_ref.data()[k*3 + 0] << "\t" << x_ref.data()[k*3 + 1] << std::endl;
486 gf_s.set_quad(k);
487 gf_f.set_quad(k);
488 // compute indeces and declare local storage
489 int eN_k = eN * nQuadraturePoints_element + k;
490 int eN_k_3d = eN_k * 3;
491 double h_phi;
492 double u = 0.0;
493 double grad_u[nSpace];
494 double ua = 0.0;
495 double grad_ua[nSpace];
496 double ub = 0.0;
497 double grad_ub[nSpace];
498 double uja = 0.0;
499 double grad_uja[nSpace];
500 double ujb = 0.0;
501 double grad_ujb[nSpace];
502 double m = 0.0;
503 double dm = 0.0;
504 double f[nSpace];
505 double df[nSpace];
506 double f_s[nSpace] = {0., 0.};
507 double df_s[nSpace] = {0., 0.};
508 double ham_s = 0.0;
509 double dham_s[nSpace] = {0., 0.};
510 double f_f[nSpace] = {0., 0.};
511 double df_f[nSpace] = {0., 0.};
512 double ham_f = 0.0;
513 double dham_f[nSpace] = {0., 0.};
514 double m_t = 0.0;
515 double dm_t = 0.0;
516 double pdeResidual_u = 0.0;
517 double Lstar_u[nDOF_test_element];
518 double subgridError_u = 0.0;
519 double tau = 0.0;
520 double tau0 = 0.0;
521 double tau1 = 0.0;
522 double numDiff0 = 0.0;
523 double numDiff1 = 0.0;
524 double *a = NULL;
525 double r = 0.0;
526 double r_s = 0.0;
527 double dr_s = 0.0;
528 double r_f = 0.0;
529 double dr_f = 0.0;
530 double jac[nSpace * nSpace];
531 double jacDet;
532 double jacInv[nSpace * nSpace];
533 double u_grad_trial[nDOF_trial_element * nSpace];
534 double u_test_dV[nDOF_trial_element];
535 double u_grad_test_dV[nDOF_test_element * nSpace];
536 double ua_grad_trial[nDOF_trial_element * nSpace];
537 double ua_test_dV[nDOF_trial_element];
538 double ua_grad_test_dV[nDOF_test_element * nSpace];
539 double ub_grad_trial[nDOF_trial_element * nSpace];
540 double ub_test_dV[nDOF_trial_element];
541 double ub_grad_test_dV[nDOF_test_element * nSpace];
542 double dV;
543 double x;
544 double y;
545 double z;
546 double G[nSpace * nSpace];
547 double G_dd_G;
548 double tr_G;
549 //
550 // compute solution and gradients at quadrature points
551 //
552 ck.calculateMapping_element(eN,
553 k,
554 mesh_dof.data(),
555 mesh_l2g.data(),
556 mesh_trial_ref.data(),
557 mesh_grad_trial_ref.data(),
558 jac,
559 jacDet,
560 jacInv,
561 x,
562 y,
563 z);
564 ck.calculateH_element(eN,
565 k,
566 nodeDiametersArray.data(),
567 mesh_l2g.data(),
568 mesh_trial_ref.data(),
569 h_phi);
570 // get the physical integration weight
571 dV = fabs(jacDet) * dV_ref.data()[k];
572 // get the metric tensor and friends
573 ck.calculateG(jacInv, G, G_dd_G, tr_G);
574 // get the trial function gradients
575 // std::cout << "Calculating gradTrialFromRef from calculateElementResidual()" << std::endl;
576 ck.gradTrialFromRef(&u_grad_trial_ref.data()[k * nDOF_trial_element * nSpace], jacInv, u_grad_trial);
577 // get the solution
578 // std::cout << "element_u data: " << std::endl;
579 ck.valFromElementDOF(element_u.data(), &u_trial_ref.data()[k * nDOF_trial_element], u);
580 // get the solution gradients
581 ck.gradFromElementDOF(element_u.data(), u_grad_trial, grad_u);
582 // precalculate test function products with integration weights
583 for (int j = 0; j < nDOF_trial_element; j++)
584 {
585 u_test_dV[j] = u_test_ref.data()[k * nDOF_trial_element + j] * dV;
586 for (int I = 0; I < nSpace; I++)
587 {
588 u_grad_test_dV[j * nSpace + I] = u_grad_trial[j * nSpace + I] * dV; // cek warning won't work for Petrov-Galerkin
589 //PGIFEM
590 //ua_grad_test_dV[j * nSpace + I] = u_grad_trial[j * nSpace + I] * dV; // cek warning won't work for Petrov-Galerkin
591 //ub_grad_test_dV[j * nSpace + I] = u_grad_trial[j * nSpace + I] * dV; // cek warning won't work for Petrov-Galerkin
592 }
593 }
594 if (icase_f == 0)
595 {
596 double va[nDOF_trial_element], va_grad_trial[nDOF_trial_element * nSpace], vb[nDOF_trial_element], vb_grad_trial[nDOF_trial_element * nSpace];
597 for (int i = 0; i < nDOF_trial_element; i++)
598 {
599 va[i] = gf_f.VA(i);
600 //assert(fabs(va[i] - u_trial_ref.data()[k * nDOF_trial_element+i])< 1.0e-8);
601 va_grad_trial[i * nSpace + 0] = gf_f.VA_x(i);
602 //assert(fabs(va_grad_trial[i * nSpace + 0] - u_grad_trial[i * nSpace + 0]) < 1.0e-8);
603 va_grad_trial[i * nSpace + 1] = gf_f.VA_y(i);
604 //assert(fabs(va_grad_trial[i * nSpace + 1] - u_grad_trial[i * nSpace + 1]) < 1.0e-8);
605 vb[i] = gf_f.VB(i);
606 //assert(fabs(vb[i] - u_trial_ref.data()[k * nDOF_trial_element+i])< 1.0e-8);
607 vb_grad_trial[i * nSpace + 0] = gf_f.VB_x(i);
608 //assert(fabs(vb_grad_trial[i * nSpace + 0] - u_grad_trial[i * nSpace + 0]) < 1.0e-8);
609 vb_grad_trial[i * nSpace + 1] = gf_f.VB_y(i);
610 //assert(fabs(vb_grad_trial[i * nSpace + 0] - u_grad_trial[i * nSpace + 0]) < 1.0e-8);
611 // std::cout << "\ni: " << i << ", va: " << va[i] << ", vb: " << vb[i] << std::endl;
612 // std::cout << "i: " << i << ", va_x: " << va_grad_trial[i * nSpace + 0] << ", va_y: " << va_grad_trial[i * nSpace + 1] << std::endl;
613 // std::cout << "i: " << i << ", vb_x: " << vb_grad_trial[i * nSpace + 0] << ", vb_y: " << vb_grad_trial[i * nSpace + 1] << std::endl;
614 }
615 // std::cout << "--------------------------------------------------------------------------------" << std::endl;
616 // std::cout << nDOF_trial_element << std::endl;
617 // std::cout << "va: " << va[0] << ", " << va[1] << ", " << va[2] << ", " << va[3] << ", " << va[4] << ", " << va[5] << std::endl;
618 // std::cout << "vb: " << vb[0] << ", " << vb[1] << ", " << vb[2] << ", " << vb[3] << ", " << vb[4] << ", " << vb[5] << std::endl;
619
620
621 //
622 //
623 // std::cout << "Calculating ua: "<< std::endl;
624 ck.valFromElementDOF(element_u.data(), va, ua);
625 ck.gradFromElementDOF(element_u.data(), va_grad_trial, grad_ua);
626 // std::cout << "Calculating ub: "<< std::endl;
627 ck.valFromElementDOF(element_u.data(), vb, ub);
628 ck.gradFromElementDOF(element_u.data(), vb_grad_trial, grad_ub);
629 // std::cout << "Calculating uja: "<< std::endl;
630 ck.valFromElementDOF(JA, va, uja);
631 ck.gradFromElementDOF(JA, va_grad_trial, grad_uja);
632 // std::cout << "Calculating ujb: "<< std::endl;
633 ck.valFromElementDOF(JB, vb, ujb);
634 ck.gradFromElementDOF(JB, vb_grad_trial, grad_ujb);
635 for (int i = 0; i < nDOF_test_element; i++)
636 {
637 ua_test_dV[i] = va[i] * dV;
638 ub_test_dV[i] = vb[i] * dV;
639 for (int I = 0; I < nSpace; I++)
640 {
641 ua_grad_test_dV[i * nSpace + I] = va_grad_trial[i * nSpace + I] * dV;
642 ub_grad_test_dV[i * nSpace + I] = vb_grad_trial[i * nSpace + I] * dV;
643 }
644 }
645 for (int j = 0; j < nDOF_trial_element; j++)
646 {
647 for (int I = 0; I < nSpace; I++)
648 {
649 ua_grad_trial[j * nSpace + I] = va_grad_trial[j * nSpace + I];
650 ub_grad_trial[j * nSpace + I] = vb_grad_trial[j * nSpace + I];
651 }
652 }
653 if (PG)
654 {
655 // Petrov-Galerkin: test with ordinary P1 hat functions instead of
656 // the enriched va/vb branches; trial/solution side is unchanged.
657 for (int i = 0; i < nDOF_test_element; i++)
658 {
659 ua_test_dV[i] = u_test_dV[i];
660 ub_test_dV[i] = u_test_dV[i];
661 for (int I = 0; I < nSpace; I++)
662 {
663 ua_grad_test_dV[i * nSpace + I] = u_grad_test_dV[i * nSpace + I];
664 ub_grad_test_dV[i * nSpace + I] = u_grad_test_dV[i * nSpace + I];
665 }
666 }
667 }
668 }
669 //
670 // calculate pde coefficients at quadrature points
671 //
672 // evaluateCoefficients();
673 // just set from pre-evaluated quadrature point values for now
674 a = &q_a.data()[eN_k * sd_rowptr.data()[nSpace]];
675 r = q_r.data()[eN_k];
676 for (int I = 0; I < nSpace; I++)
677 {
678 f[I] = q_v.data()[eN_k * nSpace + I] * u;
679 df[I] = q_v.data()[eN_k * nSpace + I];
680 }
681 const double H_s = gf_s.H(0., 0.);
682 const double D_s = gf_s.D(0., 0.);
683 if (embeddedBoundary)
684 {
685 double level_set_normal[nSpace];
686 double sign = 0.0;
687 double norm_exact = 0.0, norm_cut = 0.0;
688 for (int I = 0; I < nSpace; I++)
689 {
690 sign += embeddedBoundary_normal_q.data()[eN_k_3d + I] * gf_s.get_normal()[I];
691 level_set_normal[I] = gf_s.get_normal()[I];
692 norm_cut += level_set_normal[I] * level_set_normal[I];
693 norm_exact += embeddedBoundary_normal_q.data()[eN_k_3d + I] * embeddedBoundary_normal_q.data()[eN_k_3d + I];
694 }
695 assert(std::fabs(1.0 - norm_cut) < 1.0e-8);
696 assert(std::fabs(1.0 - norm_exact) < 1.0e-8);
697 if (sign < 0.0)
698 for (int I = 0; I < nSpace; I++)
699 level_set_normal[I] *= -1.0;
700 updateEmbeddedBoundaryTerms(embeddedBoundary_penalty / h_phi, // penalty,
701 dV,
702 level_set_normal,
703 embeddedBoundary_u_q.data()[eN_k],
704 u,
705 grad_u,
706 a[0], // assume scalar diffusion for now
707 r_s,
708 dr_s,
709 ham_s,
710 dham_s,
711 f_s,
712 df_s,
713 D_s);
714 }
715 const double ImH_f = gf_f.ImH(0., 0.);
716 const double H_f = gf_f.H(0., 0.);
717 const double D_f = gf_f.D(0., 0.);
718 // if ( H_s*ImH_f != 0.0 || D_s != 0.0 || D_f != 0.0) //for two embedded interfaces
719 if (H_s != 0.0 || D_s != 0.0 || D_f != 0.0) // for one embedded interface and one immersed interface
720 {
721 element_active = true;
722 elementIsActive[eN] = true;
723 }
724 if (immersedBoundary)
725 {
726 double level_set_normal[nSpace];
727 double sign = 0.0;
728 double norm_exact = 0.0, norm_cut = 0.0;
729 for (int I = 0; I < nSpace; I++)
730 {
731 sign += immersedBoundary_normal_q.data()[eN_k_3d + I] * gf_f.get_normal()[I];
732 level_set_normal[I] = gf_f.get_normal()[I];
733 norm_cut += level_set_normal[I] * level_set_normal[I];
734 norm_exact += immersedBoundary_normal_q.data()[eN_k_3d + I] * immersedBoundary_normal_q.data()[eN_k_3d + I];
735 }
736 assert(std::fabs(1.0 - norm_cut) < 1.0e-8);
737 assert(std::fabs(1.0 - norm_exact) < 1.0e-8);
738 if (sign < 0.0)
739 for (int I = 0; I < nSpace; I++)
740 level_set_normal[I] *= -1.0;
741 updateImmersedBoundaryTerms(immersedBoundary_penalty / h_phi, // penalty,
742 dV,
743 level_set_normal,
744 x,
745 y,
746 z,
747 immersedBoundary_u_q.data()[eN_k],
748 u,
749 grad_u,
750 a[0], // assume scalar diffusion for now
751 r_f,
752 dr_f,
753 ham_f,
754 dham_f,
755 f_f,
756 df_f,
757 immersedBoundary_fluxJump_q.data()[eN_k],
758 &immersedBoundary_fluxJumpVector_q.data()[eN_k_3d],
759 D_f);
760 }
761 //
762 // moving mesh
763 //
764 /* double mesh_velocity[3]; */
765 /* mesh_velocity[0] = xt; */
766 /* mesh_velocity[1] = yt; */
767 /* mesh_velocity[2] = zt; */
768 /* for (int I=0;I<nSpace;I++) */
769 /* { */
770 /* f[I] -= MOVING_DOMAIN*m*mesh_velocity[I]; */
771 /* df[I] -= MOVING_DOMAIN*dm*mesh_velocity[I]; */
772 /* } */
773 //
774 // calculate time derivative at quadrature points
775 //
776 /* ck.bdf(alphaBDF, */
777 /* q_m_betaBDF.data()[eN_k], */
778 /* m, */
779 /* dm, */
780 /* m_t, */
781 /* dm_t); */
782 //
783 // calculate subgrid error (strong residual and adjoint)
784 //
785 // calculate strong residual
786 pdeResidual_u = ck.Advection_strong(df, grad_u) + ck.Reaction_strong(r); // ck.Mass_strong(m_t) + ck.Advection_strong(df,grad_u) + ck.Reaction_strong(r);
787 // calculate adjoint
788 for (int i = 0; i < nDOF_test_element; i++)
789 {
790 // int eN_k_i_nSpace = (eN_k*nDOF_trial_element+i)*nSpace;
791 // Lstar_u[i] = ck.Advection_adjoint(df,&u_grad_test_dV.data()[eN_k_i_nSpace]);
792 int i_nSpace = i * nSpace;
793 Lstar_u[i] = ck.Advection_adjoint(df, &u_grad_test_dV[i_nSpace]);
794 }
795 // calculate tau and tau*Res
796 calculateSubgridError_tau(elementDiameter.data()[eN], dm_t, df, cfl.data()[eN_k], tau0);
797 calculateSubgridError_tau(Ct_sge,
798 G,
799 dm_t,
800 df,
801 tau1,
802 cfl.data()[eN_k]);
803
804 tau = useMetrics * tau1 + (1.0 - useMetrics) * tau0;
805
806 subgridError_u = -tau * pdeResidual_u;
807 //
808 // calculate shock capturing diffusion
809 //
810 ck.calculateNumericalDiffusion(shockCapturingDiffusion, elementDiameter.data()[eN], pdeResidual_u, grad_u, numDiff0);
811 ck.calculateNumericalDiffusion(shockCapturingDiffusion, sc_uref, sc_alpha, G, G_dd_G, pdeResidual_u, grad_u, numDiff1);
812 q_numDiff_u.data()[eN_k] = useMetrics * numDiff1 + (1.0 - useMetrics) * numDiff0;
813 //
814 // update element residual
815 //
816 // Leveque & Li 1994, Examples 1, 3, 4, PWC, PWL, PWQ, PWcubic
817 double a_loc[nSpace * nSpace];
818 for (int I = 0; I < nSpace * nSpace; I++) a_loc[I] = 0.0;
819
820 for (int i = 0; i < nDOF_test_element; i++)
821 {
822 int i_nSpace = i * nSpace;
823 if (icase_f == 0)
824 {
825 if(!gf_f.exact.edge && !gf_f.exact.corner)//full cut or cut on boundary of negative cell
826 {
827 for (int I = 0; I < nSpace; I++) a_loc[I * nSpace + I] = mua;
828 elementResidual_u.data()[i] += ImH_f * H_s * (ck.Advection_weak(f, &ua_grad_test_dV[i_nSpace]) +
829 ck.Diffusion_weak(sd_rowptr.data(), sd_colind.data(), a_loc, grad_ua, &ua_grad_test_dV[i_nSpace]) +
830 ck.Diffusion_weak(sd_rowptr.data(), sd_colind.data(), a_loc, grad_uja, &ua_grad_test_dV[i_nSpace]) +
831 ck.Reaction_weak(r, ua_test_dV[i]) +
832 ck.NumericalDiffusion(q_numDiff_u_last.data()[eN_k], grad_ua, &ua_grad_test_dV[i_nSpace]));
833
834 for (int I = 0; I < nSpace; I++) a_loc[I * nSpace + I] = mub;
835 elementResidual_u.data()[i] += H_f * H_s * (ck.Advection_weak(f, &ub_grad_test_dV[i_nSpace]) +
836 ck.Diffusion_weak(sd_rowptr.data(), sd_colind.data(), a_loc, grad_ub, &ub_grad_test_dV[i_nSpace]) +
837 ck.Diffusion_weak(sd_rowptr.data(), sd_colind.data(), a_loc, grad_ujb, &ub_grad_test_dV[i_nSpace]) +
838 ck.Reaction_weak(r, ub_test_dV[i]) +
839 ck.NumericalDiffusion(q_numDiff_u_last.data()[eN_k], grad_ub, &ub_grad_test_dV[i_nSpace]));
840 }
841 else if (gf_f.exact.edge == -1 || gf_f.exact.corner == -1)
842 {
843 for (int I = 0; I < nSpace; I++) a_loc[I * nSpace + I] = mua;
844 elementResidual_u.data()[i] += ImH_f * H_s * (ck.Advection_weak(f, &ua_grad_test_dV[i_nSpace]) +
845 ck.Diffusion_weak(sd_rowptr.data(), sd_colind.data(), a_loc, grad_ua, &ua_grad_test_dV[i_nSpace]) +
846 ck.Diffusion_weak(sd_rowptr.data(), sd_colind.data(), a_loc, grad_uja, &ua_grad_test_dV[i_nSpace]) +
847 ck.Reaction_weak(r, ua_test_dV[i]) +
848 ck.NumericalDiffusion(q_numDiff_u_last.data()[eN_k], grad_ua, &ua_grad_test_dV[i_nSpace]));
849 }
850 else if (gf_f.exact.edge == 1 || gf_f.exact.corner == 1)
851 {
852 for (int I = 0; I < nSpace; I++) a_loc[I * nSpace + I] = mub;
853 elementResidual_u.data()[i] += H_f * H_s * (ck.Advection_weak(f, &ub_grad_test_dV[i_nSpace]) +
854 ck.Diffusion_weak(sd_rowptr.data(), sd_colind.data(), a_loc, grad_ub, &ub_grad_test_dV[i_nSpace]) +
855 ck.Diffusion_weak(sd_rowptr.data(), sd_colind.data(), a_loc, grad_ujb, &ub_grad_test_dV[i_nSpace]) +
856 ck.Reaction_weak(r, ub_test_dV[i]) +
857 ck.NumericalDiffusion(q_numDiff_u_last.data()[eN_k], grad_ub, &ub_grad_test_dV[i_nSpace]));
858 }
859 else assert(false && "Invalid gf_f.exact.edge/corner values. Should be -1, 0 or +1.");
860 }
861 else
862 {
863
864 elementResidual_u.data()[i] += H_s * (ck.Advection_weak(f, &u_grad_test_dV[i_nSpace]) +
865 ck.Diffusion_weak(sd_rowptr.data(), sd_colind.data(), a, grad_u, &u_grad_test_dV[i_nSpace]) +
866 ck.Reaction_weak(r, u_test_dV[i]) +
867 ck.SubgridError(subgridError_u, Lstar_u[i]) +
868 ck.NumericalDiffusion(q_numDiff_u_last.data()[eN_k], grad_u, &u_grad_test_dV[i_nSpace]));
869 }
870 if (embeddedBoundary)
871 {
872 if (gf_s.exact.edge >= 0 && !gf_s.exact.corner)
873 {
874 elementResidual_u.data()[i] += (ck.Advection_weak(f_s, &u_grad_test_dV[i_nSpace]) +
875 ck.Reaction_weak(r_s, u_test_dV[i]) +
876 ck.Hamiltonian_weak(ham_s, u_test_dV[i]));
877 }
878 }
879 if (immersedBoundary)
880 {
881 if (gf_f.exact.edge >= 0 && !gf_f.exact.corner)
882 {
883 elementResidual_u.data()[i] += (ck.Advection_weak(f_f, &u_grad_test_dV[i_nSpace]) +
884 ck.Reaction_weak(r_f, u_test_dV[i]) +
885 ck.Hamiltonian_weak(ham_f, u_test_dV[i]));
886 }
887 }
888 }
889 double L2_contrib = 0.0;
890 if (icase_f == 0)
891 {
892 double sol_in = q_u_exact_inner.data()[eN_k];
893 double err_in = fabs(ua + uja - sol_in);
894 L2_contrib += ImH_f * err_in * err_in * dV;
895 double sol_out = q_u_exact_outer.data()[eN_k];
896 double err_out = fabs(ub + ujb - sol_out);
897 L2_contrib += H_f * err_out * err_out * dV;
898 if (ImH_f >= H_f)
899 Linfty_error = std::max(Linfty_error, err_in);
900 else
901 Linfty_error = std::max(Linfty_error, err_out);
902 }
903 else
904 {
905 if (icase_f == -1)
906 {
907 double sol = q_u_exact_inner.data()[eN_k];
908 double err = fabs(u - sol);
909 L2_contrib += err * err * dV;
910 Linfty_error = std::max(Linfty_error, err);
911 }
912 if (icase_f == 1)
913 {
914 double sol = q_u_exact_outer.data()[eN_k];
915 double err = fabs(u - sol);
916 L2_contrib += err * err * dV;
917 Linfty_error = std::max(Linfty_error, err);
918 }
919 }
920 L2_error += L2_contrib;
921 }
922 }
923
924 void calculateResidual(arguments_dict &args)
925 {
926 xt::pyarray<double> &mesh_trial_ref = args.array<double>("mesh_trial_ref");
927 xt::pyarray<double> &mesh_grad_trial_ref = args.array<double>("mesh_grad_trial_ref");
928 xt::pyarray<double> &mesh_dof = args.array<double>("mesh_dof");
929 xt::pyarray<int> &mesh_l2g = args.array<int>("mesh_l2g");
930 xt::pyarray<double> &dV_ref = args.array<double>("dV_ref");
931 xt::pyarray<double> &u_trial_ref = args.array<double>("u_trial_ref");
932 xt::pyarray<double> &u_grad_trial_ref = args.array<double>("u_grad_trial_ref");
933 xt::pyarray<double> &u_test_ref = args.array<double>("u_test_ref");
934 xt::pyarray<double> &u_grad_test_ref = args.array<double>("u_grad_test_ref");
935 xt::pyarray<double> &elementDiameter = args.array<double>("elementDiameter");
936 xt::pyarray<double> &cfl = args.array<double>("cfl");
937 double Ct_sge = args.scalar<double>("Ct_sge");
938 double sc_uref = args.scalar<double>("sc_uref");
939 double sc_alpha = args.scalar<double>("sc_alpha");
940 double useMetrics = args.scalar<double>("useMetrics");
941 xt::pyarray<double> &mesh_trial_trace_ref = args.array<double>("mesh_trial_trace_ref");
942 xt::pyarray<double> &mesh_grad_trial_trace_ref = args.array<double>("mesh_grad_trial_trace_ref");
943 xt::pyarray<double> &dS_ref = args.array<double>("dS_ref");
944 xt::pyarray<double> &u_trial_trace_ref = args.array<double>("u_trial_trace_ref");
945 xt::pyarray<double> &u_grad_trial_trace_ref = args.array<double>("u_grad_trial_trace_ref");
946 xt::pyarray<double> &u_test_trace_ref = args.array<double>("u_test_trace_ref");
947 xt::pyarray<double> &u_grad_test_trace_ref = args.array<double>("u_grad_test_trace_ref");
948 xt::pyarray<double> &normal_ref = args.array<double>("normal_ref");
949 xt::pyarray<double> &boundaryJac_ref = args.array<double>("boundaryJac_ref");
950 int nElements_global = args.scalar<int>("nElements_global");
951 xt::pyarray<int> &u_l2g = args.array<int>("u_l2g");
952 xt::pyarray<double> &u_dof = args.array<double>("u_dof");
953 xt::pyarray<int> &sd_rowptr = args.array<int>("sd_rowptr");
954 xt::pyarray<int> &sd_colind = args.array<int>("sd_colind");
955 xt::pyarray<double> &q_a = args.array<double>("q_a");
956 xt::pyarray<double> &q_v = args.array<double>("q_v");
957 xt::pyarray<double> &q_r = args.array<double>("q_r");
958 int lag_shockCapturing = args.scalar<int>("lag_shockCapturing");
959 double shockCapturingDiffusion = args.scalar<double>("shockCapturingDiffusion");
960 xt::pyarray<double> &q_numDiff_u = args.array<double>("q_numDiff_u");
961 xt::pyarray<double> &q_numDiff_u_last = args.array<double>("q_numDiff_u_last");
962 int offset_u = args.scalar<int>("offset_u");
963 int stride_u = args.scalar<int>("stride_u");
964 xt::pyarray<double> &globalResidual = args.array<double>("globalResidual");
965 int nExteriorElementBoundaries_global = args.scalar<int>("nExteriorElementBoundaries_global");
966 xt::pyarray<int> &exteriorElementBoundariesArray = args.array<int>("exteriorElementBoundariesArray");
967 xt::pyarray<int> &elementBoundaryElementsArray = args.array<int>("elementBoundaryElementsArray");
968 xt::pyarray<int> &elementBoundaryLocalElementBoundariesArray = args.array<int>("elementBoundaryLocalElementBoundariesArray");
969 xt::pyarray<double> &ebqe_a = args.array<double>("ebqe_a");
970 xt::pyarray<double> &ebqe_v = args.array<double>("ebqe_v");
971 xt::pyarray<int> &isDOFBoundary_u = args.array<int>("isDOFBoundary_u");
972 xt::pyarray<double> &ebqe_bc_u_ext = args.array<double>("ebqe_bc_u_ext");
973 xt::pyarray<int> &isDiffusiveFluxBoundary_u = args.array<int>("isDiffusiveFluxBoundary_u");
974 xt::pyarray<int> &isAdvectiveFluxBoundary_u = args.array<int>("isAdvectiveFluxBoundary_u");
975 xt::pyarray<double> &ebqe_bc_flux_u_ext = args.array<double>("ebqe_bc_flux_u_ext");
976 xt::pyarray<double> &ebqe_bc_advectiveFlux_u_ext = args.array<double>("ebqe_bc_advectiveFlux_u_ext");
977 xt::pyarray<double> &ebqe_penalty_ext = args.array<double>("ebqe_penalty_ext");
978 const bool embeddedBoundary = args.scalar<int>("embeddedBoundary");
979 const double embeddedBoundary_penalty = args.scalar<double>("embeddedBoundary_penalty");
980 const double embeddedBoundary_ghost_penalty = args.scalar<double>("embeddedBoundary_ghost_penalty");
981 xt::pyarray<double> &embeddedBoundary_sdf_nodes = args.array<double>("embeddedBoundary_sdf_nodes");
982 xt::pyarray<double> &embeddedBoundary_sdf_q = args.array<double>("embeddedBoundary_sdf_q");
983 xt::pyarray<double> &embeddedBoundary_normal_q = args.array<double>("embeddedBoundary_normal_q");
984 xt::pyarray<double> &embeddedBoundary_u_q = args.array<double>("embeddedBoundary_u_q");
985 const bool immersedBoundary = args.scalar<int>("immersedBoundary");
986 const double immersedBoundary_penalty = args.scalar<double>("immersedBoundary_penalty");
987 const double immersedSCIFEM_switch = args.scalar<double>("immersedSCIFEM_switch");
988 const double immersedSCIFEM_penalty = args.scalar<double>("immersedSCIFEM_penalty");
989 const bool PG = args.scalar<int>("PG");
990 xt::pyarray<double> &immersedBoundary_sdf_nodes = args.array<double>("immersedBoundary_sdf_nodes");
991 xt::pyarray<double> &immersedBoundary_sdf_q = args.array<double>("immersedBoundary_sdf_q");
992 xt::pyarray<double> &immersedBoundary_normal_q = args.array<double>("immersedBoundary_normal_q");
993 xt::pyarray<double> &immersedBoundary_u_q = args.array<double>("immersedBoundary_u_q");
994 xt::pyarray<double> &immersedBoundary_fluxJump_q = args.array<double>("immersedBoundary_fluxJump_q");
995 xt::pyarray<double> &immersedBoundary_fluxJumpVector_q = args.array<double>("immersedBoundary_fluxJumpVector_q");
996 xt::pyarray<double> &immersedBoundary_solutionJump_nodes = args.array<double>("immersedBoundary_solutionJump_nodes");
997 xt::pyarray<double> &isActiveDOF = args.array<double>("isActiveDOF");
998 const double eb_adjoint_sigma = args.scalar<double>("eb_adjoint_sigma");
999 xt::pyarray<double> &x_ref = args.array<double>("x_ref");
1000 xt::pyarray<double> &xB_ref = args.array<double>("xB_ref");
1001 xt::pyarray<int> &elementBoundariesArray = args.array<int>("elementBoundariesArray");
1002 const int nElementBoundaries_owned = args.scalar<int>("nElementBoundaries_owned");
1003 xt::pyarray<double> &elementBoundaryDiameter = args.array<double>("elementBoundaryDiameter");
1004 xt::pyarray<double> &nodeDiametersArray = args.array<double>("nodeDiametersArray");
1005 xt::pyarray<double> &L2_error = args.array<double>("L2_error");
1006 xt::pyarray<double> &Linfty_error = args.array<double>("Linfty_error");
1007 const double mua = args.scalar<double>("mua");
1008 const double mub = args.scalar<double>("mub");
1009 const double jf = args.scalar<double>("jf");
1010 xt::pyarray<double> &q_u_exact_inner = args.array<double>("q_u_exact_inner");
1011 xt::pyarray<double> &q_u_exact_outer = args.array<double>("q_u_exact_outer");
1012 const bool recomputeIFEMGeometry = args.scalar<int>("recomputeIFEMGeometry");
1013 ensureIFEMCacheSized(nElements_global); // also (re)asserts useExact = true on (re)allocation
1014 if (recomputeIFEMGeometry)
1015 ifemGeometryGeneration++;
1016 ifem_boundaries.clear();
1017 ifem_boundary_elements.clear();
1018 cutfem_boundaries.clear();
1019 cutfem_boundary_elements.clear();
1020 elementIsActive.resize(nElements_global);
1021
1022 //
1023 // loop over elements to compute volume integrals and load them into element and global residual
1024 //
1025 // eN is the element index
1026 // eN_k is the quadrature point index for a scalar
1027 // eN_k_nSpace is the quadrature point index for a vector
1028 // eN_i is the element test function index
1029 // eN_j is the element trial function index
1030 // eN_k_j is the quadrature point index for a trial function
1031 // eN_k_i is the quadrature point index for a trial function
1032 for (int eN = 0; eN < nElements_global; eN++)
1033 {
1034 // std::cout << "########################\n element: " << eN << " \n########################" << std::endl;
1035 // declare local storage for element residual and initialize
1036 // double elementResidual_u[nDOF_test_element],element_u[nDOF_trial_element];
1037 auto elementResidual_u = xt::pyarray<double>::from_shape({nDOF_test_element});
1038 auto element_u = xt::pyarray<double>::from_shape({nDOF_trial_element});
1039 bool element_active = false;
1040 elementIsActive[eN] = false;
1041 for (int i = 0; i < nDOF_trial_element; i++)
1042 {
1043 int eN_i = eN * nDOF_trial_element + i;
1044 element_u.data()[i] = u_dof.data()[u_l2g.data()[eN_i]];
1045 // std::cout << "element_u[" << i << "]:" << element_u.data()[i] << std::endl;
1046 } // i
1047 double element_phi_s[nDOF_trial_element];
1048 for (int j = 0; j < nDOF_trial_element; j++)
1049 {
1050 int eN_j = eN * nDOF_trial_element + j;
1051 element_phi_s[j] = embeddedBoundary_sdf_nodes.data()[u_l2g.data()[eN_j]];
1052 }
1053 double element_phi_f[nDOF_trial_element];
1054 for (int j = 0; j < nDOF_trial_element; j++)
1055 {
1056 int eN_j = eN * nDOF_trial_element + j;
1057 element_phi_f[j] = immersedBoundary_sdf_nodes.data()[u_l2g.data()[eN_j]];
1058 }
1059 // std::cout << std::endl;
1060 double element_nodes[nDOF_trial_element * 3];
1061 for (int i = 0; i < nDOF_trial_element; i++)
1062 {
1063 int eN_i = eN * nDOF_trial_element + i;
1064 for (int I = 0; I < 3; I++)
1065 // element_nodes[i * 3 + I] = mesh_dof.data()[mesh_l2g.data()[eN_i] * 3 + I];
1066 element_nodes[i * 3 + I] = mesh_dof.data()[u_l2g.data()[eN_i] * 3 + I];
1067 // std::cout << "element node[" << i << "]:" << element_nodes[i * 3 + 0] << " " << element_nodes[i * 3 + 1] << " " << element_nodes[i * 3 + 2] << std::endl;
1068 // std::cout << "element node[" << i << "]:" << element_nodes[i * 3 + 0] << " " << element_nodes[i * 3 + 1] << " " << element_nodes[i * 3 + 2] << std::endl;
1069 // std::cout << "element phi_f[" << i << "]:" << element_phi_f[i] << std::endl << std::endl;
1070 } // i
1071 if (gf_s_interior_gen[eN] != ifemGeometryGeneration)
1072 {
1073 gf_s_interior_icase[eN] = gf_s_cache[eN].calculate(element_phi_s, element_nodes, x_ref.data(), false);
1074 gf_s_interior_gen[eN] = ifemGeometryGeneration;
1075 }
1076 int icase_s = gf_s_interior_icase[eN];
1077 if (icase_s == 0)
1078 {
1079 // only works for simplices
1080 for (int ebN_element = 0; ebN_element < nDOF_mesh_trial_element; ebN_element++)
1081 {
1082 const int ebN = elementBoundariesArray.data()[eN * nDOF_mesh_trial_element + ebN_element];
1083 // internal and actually a cut edge
1084 if (elementBoundaryElementsArray.data()[ebN * 2 + 1] != -1 && (ebN < nElementBoundaries_owned))
1085 cutfem_boundaries.insert(ebN);
1086 }
1087 }
1088 if (gf_f_interior_gen[eN] != ifemGeometryGeneration)
1089 {
1090 gf_f_interior_icase[eN] = gf_f_cache[eN].calculate(element_phi_f, element_nodes, x_ref.data(), mua, mub, jf, false, false);
1091 gf_f_interior_gen[eN] = ifemGeometryGeneration;
1092 }
1093 int icase_f = gf_f_interior_icase[eN];
1094 double JA[nDOF_trial_element];
1095 double JB[nDOF_trial_element];
1096 std::fill(JA, JA + nDOF_trial_element, 0.0);
1097 std::fill(JB, JB + nDOF_trial_element, 0.0);
1098 if (icase_f == 0)
1099 {
1100 // std::cout << "Active element" << std::endl;
1101 // only works for simplices
1102 for (int ebN_element = 0; ebN_element < nDOF_mesh_trial_element; ebN_element++)
1103 {
1104 // std::cout << "ebN_element =" << ebN_element << "\t ebN=" << ebN << std::endl;
1105 // internal and actually a cut edge
1106 const int ebN = elementBoundariesArray.data()[eN * nDOF_mesh_trial_element + ebN_element];
1107 // if (elementBoundaryElementsArray.data()[ebN * 2 + 1] != -1 && (ebN < nElementBoundaries_owned)) // This gives all the internal edges instead of just the cut edges.
1108
1109 // indexing convention for P1: edge opposite to node i is given by (i+1)%3 and (i+2)%3 nodes.
1110 // nathawani: P2 needs different indexing convention.
1111 // Do we need corner cases? (<=0) or just (<0) for cut edge detection?
1112 // The SCIFEM facet term is an inter-element consistency term: it needs the
1113 // trace from BOTH elements sharing the face. On an exterior face the second
1114 // entry of elementBoundaryElementsArray is -1, and the loops below index
1115 // elementIsActive[] and gf_f_cache[] with it before any validity test, which
1116 // is out of bounds. test=8's circle lies wholly inside the domain so this
1117 // never triggered; test=12's diagonal interface exits through the boundary
1118 // and segfaulted on arm64. Only interior faces belong in this set.
1119 if (elementBoundaryElementsArray.data()[ebN * 2 + 1] != -1 &&
1120 element_phi_f[(ebN_element + 1) % 3] * element_phi_f[(ebN_element + 2) % 3] < 0.0)
1121 {
1122 // This should give just the cut edges for simplices.
1123 ifem_boundaries.insert(ebN);
1124 }
1125 }
1126 // Leveque & Li 1994, Example 1, 2, 3, 4; Ji et. al. 2014
1127 // The jump function is, in general, not constant over the cut element, so it
1128 // must be evaluated at each node's own coordinates (not once at the cut
1129 // barycenter) -- otherwise JA/JB are constant over the element and their
1130 // interpolated gradients (grad_uja, grad_ujb) vanish identically.
1131 const double eps_phi = 1.0e-12;
1132 auto assign_jump_side = [&](int i, bool isOuter, double jump) {
1133 if (isOuter)
1134 {
1135 JA[i] = -jump;
1136 JB[i] = 0.0;
1137 }
1138 else
1139 {
1140 JA[i] = 0.0;
1141 JB[i] = jump;
1142 }
1143 };
1144 for (int i = 0; i < nDOF_trial_element; i++)
1145 {
1146
1147 int eN_i = eN * nDOF_trial_element + i;
1148 const double jump_i = immersedBoundary_solutionJump_nodes.data()[u_l2g.data()[eN_i]];
1149 if (element_phi_f[i] > 0.0)
1150 {
1151 assign_jump_side(i, true, jump_i);
1152 }
1153 else if (element_phi_f[i] <= 0.0)
1154 {
1155 assign_jump_side(i, false, jump_i);
1156 }
1157 }
1158 }
1159 else if (icase_f == -1)
1160 {
1161 }
1162 else if (icase_f == 1)
1163 {
1164 }
1165 calculateElementResidual(icase_f,
1166 mesh_trial_ref,
1167 mesh_grad_trial_ref,
1168 mesh_dof,
1169 mesh_l2g,
1170 x_ref,
1171 dV_ref,
1172 u_trial_ref,
1173 u_grad_trial_ref,
1174 u_test_ref,
1175 u_grad_test_ref,
1176 elementDiameter,
1177 elementBoundaryDiameter,
1178 nodeDiametersArray,
1179 cfl,
1180 Ct_sge,
1181 sc_uref,
1182 sc_alpha,
1183 useMetrics,
1184 mesh_trial_trace_ref,
1185 mesh_grad_trial_trace_ref,
1186 dS_ref,
1187 u_trial_trace_ref,
1188 u_grad_trial_trace_ref,
1189 u_test_trace_ref,
1190 u_grad_test_trace_ref,
1191 normal_ref,
1192 boundaryJac_ref,
1193 nElements_global,
1194 nElementBoundaries_owned,
1195 u_l2g,
1196 u_dof,
1197 sd_rowptr,
1198 sd_colind,
1199 q_a,
1200 q_v,
1201 q_r,
1202 lag_shockCapturing,
1203 shockCapturingDiffusion,
1204 q_numDiff_u,
1205 q_numDiff_u_last,
1206 offset_u, stride_u,
1207 elementResidual_u,
1208 nExteriorElementBoundaries_global,
1209 exteriorElementBoundariesArray,
1210 elementBoundariesArray,
1211 elementBoundaryElementsArray,
1212 elementBoundaryLocalElementBoundariesArray,
1213 element_u,
1214 eN,
1215 embeddedBoundary,
1216 embeddedBoundary_penalty,
1217 embeddedBoundary_normal_q,
1218 embeddedBoundary_u_q,
1219 immersedBoundary,
1220 immersedBoundary_penalty,
1221 immersedBoundary_sdf_q,
1222 immersedBoundary_normal_q,
1223 immersedBoundary_u_q,
1224 immersedBoundary_fluxJump_q,
1225 immersedBoundary_fluxJumpVector_q,
1226 immersedBoundary_solutionJump_nodes,
1227 element_phi_f,
1228 element_active,
1229 elementIsActive,
1230 JA,
1231 JB,
1232 L2_error.data()[0],
1233 Linfty_error.data()[0],
1234 mua,
1235 mub,
1236 q_u_exact_inner,
1237 q_u_exact_outer,
1238 PG);
1239 //
1240 // load element into global residual and save element residual
1241 //
1242 for (int i = 0; i < nDOF_test_element; i++)
1243 {
1244 int eN_i = eN * nDOF_test_element + i;
1245 globalResidual.data()[offset_u + stride_u * u_l2g.data()[eN_i]] += elementResidual_u.data()[i];
1246 if (element_active)
1247 isActiveDOF.data()[offset_u + stride_u * u_l2g.data()[eN_i]] = 1.0;
1248 // std::cout << "globalResidual[" << offset_u + stride_u * u_l2g.data()[eN_i] << "] += " << elementResidual_u.data()[i] << std::endl;
1249 } // i
1250 } // elements
1251 for (std::set<int>::iterator it = cutfem_boundaries.begin(); it != cutfem_boundaries.end();)
1252 {
1253 if (elementIsActive[elementBoundaryElementsArray[(*it) * 2 + 0]] && elementIsActive[elementBoundaryElementsArray[(*it) * 2 + 1]])
1254 {
1255 std::map<int, double> Dwp_Dn_jump, Dw_Dn_jump;
1256 double gamma_cutfem = embeddedBoundary_ghost_penalty, h_cutfem = elementBoundaryDiameter.data()[*it];
1257 for (int kb = 0; kb < nQuadraturePoints_elementBoundary; kb++)
1258 {
1259 double Du_Dn_jump = 0.0, dS;
1260 for (int eN_side = 0; eN_side < 2; eN_side++)
1261 {
1262 int ebN = *it,
1263 eN = elementBoundaryElementsArray.data()[ebN * 2 + eN_side];
1264 for (int i = 0; i < nDOF_test_element; i++)
1265 {
1266 Dw_Dn_jump[u_l2g.data()[eN * nDOF_test_element + i]] = 0.0;
1267 }
1268 }
1269 for (int eN_side = 0; eN_side < 2; eN_side++)
1270 {
1271 int ebN = *it,
1272 eN = elementBoundaryElementsArray.data()[ebN * 2 + eN_side],
1273 ebN_local = elementBoundaryLocalElementBoundariesArray.data()[ebN * 2 + eN_side],
1274 eN_nDOF_trial_element = eN * nDOF_trial_element,
1275 ebN_local_kb = ebN_local * nQuadraturePoints_elementBoundary + kb,
1276 ebN_local_kb_nSpace = ebN_local_kb * nSpace;
1277 double u_int = 0.0,
1278 grad_u_int[nSpace] = {0., 0.},
1279 jac_int[nSpace * nSpace],
1280 jacDet_int,
1281 jacInv_int[nSpace * nSpace],
1282 boundaryJac[nSpace * (nSpace - 1)],
1283 metricTensor[(nSpace - 1) * (nSpace - 1)],
1284 metricTensorDetSqrt,
1285 u_test_dS[nDOF_test_element],
1286 u_grad_trial_trace[nDOF_trial_element * nSpace],
1287 u_grad_test_dS[nDOF_trial_element * nSpace],
1288 normal[2], x_int, y_int, z_int, xt_int, yt_int, zt_int, integralScaling,
1289 G[nSpace * nSpace], G_dd_G, tr_G, h_phi, h_penalty, penalty,
1290 force_x, force_y, force_z, force_p_x, force_p_y, force_p_z, force_v_x, force_v_y, force_v_z, r_x, r_y, r_z;
1291 // compute information about mapping from reference element to physical element
1292 ck.calculateMapping_elementBoundary(eN,
1293 ebN_local,
1294 kb,
1295 ebN_local_kb,
1296 mesh_dof.data(),
1297 mesh_l2g.data(),
1298 mesh_trial_trace_ref.data(),
1299 mesh_grad_trial_trace_ref.data(),
1300 boundaryJac_ref.data(),
1301 jac_int,
1302 jacDet_int,
1303 jacInv_int,
1304 boundaryJac,
1305 metricTensor,
1306 metricTensorDetSqrt,
1307 normal_ref.data(),
1308 normal,
1309 x_int, y_int, z_int);
1310 dS = metricTensorDetSqrt * dS_ref.data()[kb];
1311 // compute shape and solution information
1312 // shape
1313 // std::cout << "Calculating gradTrialFromRef from calculateResidual() 1" << std::endl;
1314 ck.gradTrialFromRef(&u_grad_trial_trace_ref.data()[ebN_local_kb_nSpace * nDOF_trial_element], jacInv_int, u_grad_trial_trace);
1315 // solution and gradients
1316 ck.valFromDOF(u_dof.data(), &u_l2g.data()[eN_nDOF_trial_element], &u_trial_trace_ref.data()[ebN_local_kb * nDOF_test_element], u_int);
1317 ck.gradFromDOF(u_dof.data(), &u_l2g.data()[eN_nDOF_trial_element], u_grad_trial_trace, grad_u_int);
1318 for (int I = 0; I < nSpace; I++)
1319 {
1320 Du_Dn_jump += grad_u_int[I] * normal[I];
1321 }
1322 for (int i = 0; i < nDOF_test_element; i++)
1323 {
1324 for (int I = 0; I < nSpace; I++)
1325 Dw_Dn_jump[u_l2g.data()[eN_nDOF_trial_element + i]] += u_grad_trial_trace[i * nSpace + I] * normal[I];
1326 }
1327 } // eN_side
1328 for (std::map<int, double>::iterator w_it = Dw_Dn_jump.begin(); w_it != Dw_Dn_jump.end(); ++w_it)
1329 {
1330 int i_global = w_it->first;
1331 double Dw_Dn_jump_i = w_it->second;
1332 globalResidual.data()[offset_u + stride_u * i_global] += gamma_cutfem * h_cutfem * Du_Dn_jump * Dw_Dn_jump_i * dS;
1333 } // i
1334 } // kb
1335 ++it;
1336 }
1337 else
1338 {
1339 it = cutfem_boundaries.erase(it);
1340 }
1341 } // cutfem element boundaries
1342 // Per-face reference-ELEMENT quadrature points, taken from proteus' own trace tables.
1343 // Simplex::calculate(..., isBoundary=true) reads its xi_r argument as reference *element*
1344 // coordinates (it forms x = node0 + Jac_0*xi_r), but xB_ref holds reference *boundary*
1345 // points (t,0,0) -- passing those directly puts every face's points on the local
1346 // node0->node1 edge, i.e. correct only for local face 2. mesh_trial_trace_ref holds the P1
1347 // mesh shape functions at (face, quadrature point) and for a triangle those barycentric
1348 // functions ARE the reference coordinates (phi_0=1-xi-eta, phi_1=xi, phi_2=eta), so this
1349 // is exactly consistent with calculateMapping_elementBoundary's own convention.
1350 double xB_ref_faces[nDOF_mesh_trial_element * nQuadraturePoints_elementBoundary * 3];
1351 for (int f = 0; f < nDOF_mesh_trial_element; f++)
1352 for (int kb = 0; kb < nQuadraturePoints_elementBoundary; kb++)
1353 {
1354 const int fkb = f * nQuadraturePoints_elementBoundary + kb;
1355 const double *phi_m = &mesh_trial_trace_ref.data()[fkb * nDOF_mesh_trial_element];
1356 xB_ref_faces[fkb * 3 + 0] = phi_m[1];
1357 xB_ref_faces[fkb * 3 + 1] = phi_m[2];
1358 xB_ref_faces[fkb * 3 + 2] = 0.0;
1359 }
1360 // SCIFEM (Ji et al. 2014, eq. 3.7/3.8) symmetric consistency + adjoint-consistency terms
1361 // on cut edges: -switch*({{beta grad u}}.n_e*[[v]] + {{beta grad v}}.n_e*[[u]]).
1362 // A genuine 1D facet term -- the edge carries its own arc-length measure -- so no Dirac
1363 // surrogate is involved. The two region-wise integrands are formed separately and blended
1364 // by the edge Heaviside fit (see the prologue below), which is what makes this exact for
1365 // P1 and P2 without splitting the quadrature.
1366 for (std::set<int>::iterator it = ifem_boundaries.begin(); it != ifem_boundaries.end();)
1367 {
1368 if (elementIsActive[elementBoundaryElementsArray[(*it) * 2 + 0]] && elementIsActive[elementBoundaryElementsArray[(*it) * 2 + 1]])
1369 {
1370 int ebN = *it;
1371 int eN_s[2], ebN_local_s[2];
1372 for (int s = 0; s < 2; s++)
1373 {
1374 eN_s[s] = elementBoundaryElementsArray.data()[ebN * 2 + s];
1375 ebN_local_s[s] = elementBoundaryLocalElementBoundariesArray.data()[ebN * 2 + s];
1376 }
1377
1378 // --- Orientation map ------------------------------------------------------
1379 // The two elements sharing an interior face do NOT necessarily parametrize it in
1380 // the same direction, so quadrature index kb can denote DIFFERENT physical points
1381 // on the two sides (measured: 4 of 14 cut edges at test=8.0/refinement=1, with
1382 // exactly the mirror offsets |1-2t_k|*L of the Gauss rule). Pairing the sides by
1383 // raw kb then differences u at two different points, so [[u]] and {{beta grad u}}
1384 // pick up u's variation ALONG the edge instead of the genuine inter-element jump.
1385 // Build an explicit map instead: kmap[s][k] is side s's index for the physical
1386 // point that side 0 calls k. Matched by position, so no assumption about the
1387 // rule being symmetric or ascending.
1388 double xq[2][nQuadraturePoints_elementBoundary][3];
1389 for (int s = 0; s < 2; s++)
1390 for (int kb = 0; kb < nQuadraturePoints_elementBoundary; kb++)
1391 {
1392 double jac_t[nSpace * nSpace], jacDet_t, jacInv_t[nSpace * nSpace],
1393 bJac_t[nSpace * (nSpace - 1)], mT_t[(nSpace - 1) * (nSpace - 1)],
1394 mTDS_t, nrm_t[2], xt_ = 0.0, yt_ = 0.0, zt_ = 0.0;
1395 ck.calculateMapping_elementBoundary(eN_s[s], ebN_local_s[s], kb,
1396 ebN_local_s[s] * nQuadraturePoints_elementBoundary + kb,
1397 mesh_dof.data(), mesh_l2g.data(),
1398 mesh_trial_trace_ref.data(), mesh_grad_trial_trace_ref.data(),
1399 boundaryJac_ref.data(), jac_t, jacDet_t, jacInv_t,
1400 bJac_t, mT_t, mTDS_t, normal_ref.data(), nrm_t,
1401 xt_, yt_, zt_);
1402 xq[s][kb][0] = xt_; xq[s][kb][1] = yt_; xq[s][kb][2] = zt_;
1403 }
1404 // edge diameter: the length scale for the pairing tolerance just below,
1405 // and for the interior-penalty scaling further down (same measure the
1406 // cutfem ghost penalty uses); note gamma/h here for a VALUE jump, vs
1407 // gamma*h there for a gradient jump.
1408 const double h_edge = elementBoundaryDiameter.data()[ebN];
1409 // Match over nSpace components only: a component the mapping does not
1410 // define contributes nothing but noise to the distance, and the paired
1411 // points agree EXACTLY in the components it does define, so the winning
1412 // distance is 0 and the runner-up is a full quadrature spacing away.
1413 // Verify that outcome rather than trust it -- a mis-pick here silently
1414 // differences u at two different points, and assert() is compiled out by
1415 // the -DNDEBUG that reaches these translation units.
1416 int kmap[2][nQuadraturePoints_elementBoundary];
1417 bool kmap_used[nQuadraturePoints_elementBoundary] = {};
1418 for (int k = 0; k < nQuadraturePoints_elementBoundary; k++)
1419 {
1420 kmap[0][k] = k;
1421 int best = -1;
1422 double bestd = 1.0e300, nextd = 1.0e300;
1423 for (int j = 0; j < nQuadraturePoints_elementBoundary; j++)
1424 {
1425 double d = 0.0;
1426 for (int I = 0; I < nSpace; I++)
1427 d += std::fabs(xq[1][j][I] - xq[0][k][I]);
1428 if (d < bestd) { nextd = bestd; bestd = d; best = j; }
1429 else if (d < nextd) { nextd = d; }
1430 }
1431 if (best < 0 || kmap_used[best] || bestd >= 1.0e-8 * h_edge || nextd <= 2.0 * bestd)
1432 {
1433 std::cerr << "ADR immersedSCIFEM: quadrature points on face " << ebN
1434 << " between elements " << eN_s[0] << " and " << eN_s[1]
1435 << " do not pair up (k=" << k << ", best=" << best
1436 << ", d=" << bestd << ", runner-up=" << nextd
1437 << ", h=" << h_edge << "); the mesh is non-conforming or the"
1438 << " boundary mapping is inconsistent." << std::endl;
1439 throw std::runtime_error("ADR immersedSCIFEM: face quadrature pairing failed");
1440 }
1441 kmap_used[best] = true;
1442 kmap[1][k] = best;
1443 }
1444
1445 // --- Edge Heaviside moment fit -------------------------------------------
1446 // The facet integrand is DISCONTINUOUS at the interface crossing theta, so a fixed
1447 // whole-edge Gauss rule cannot integrate it directly. Weight the two one-sided
1448 // integrands by an edge-restricted moment-fit Heaviside instead: exact for any
1449 // polynomial integrand of degree <= nP_edge, hence exact for P1 (degree <= 1) and
1450 // P2 (degree <= 3) alike -- no linearity assumption, unlike a split-quadrature
1451 // scheme. ONE fit per edge, in side 0's parametrization: the two sides' fits are
1452 // mirror images, moment-equivalent but not pointwise equal, so they must not be
1453 // mixed. t below is therefore always measured along side 0's edge direction.
1454 double edge_nodes[2][3], phi_edge[2];
1455 {
1456 const int n0l = (ebN_local_s[0] + 1) % 3, n1l = (ebN_local_s[0] + 2) % 3;
1457 const int ln[2] = {n0l, n1l};
1458 for (int e = 0; e < 2; e++)
1459 {
1460 int eN_i = eN_s[0] * nDOF_trial_element + ln[e];
1461 phi_edge[e] = immersedBoundary_sdf_nodes.data()[u_l2g.data()[eN_i]];
1462 for (int I = 0; I < 3; I++)
1463 edge_nodes[e][I] = mesh_dof.data()[u_l2g.data()[eN_i] * 3 + I];
1464 }
1465 }
1466 double C_H_edge[nP_edge + 1];
1467 equivalent_polynomials::calculate_edge_H<nP_edge>(phi_edge[0], phi_edge[1], C_H_edge);
1468 double edge_vec[3];
1469 double edge_len2 = 0.0;
1470 for (int I = 0; I < 3; I++)
1471 {
1472 edge_vec[I] = edge_nodes[1][I] - edge_nodes[0][I];
1473 edge_len2 += edge_vec[I] * edge_vec[I];
1474 }
1475
1476 for (int k = 0; k < nQuadraturePoints_elementBoundary; k++)
1477 {
1478 double dS = 0.0;
1479 double ua_s[2], ub_s[2], flux_a_s[2], flux_b_s[2];
1480 std::map<int, double> va_m[2], vb_m[2], fta_m[2], ftb_m[2];
1481 for (int eN_side = 0; eN_side < 2; eN_side++)
1482 {
1483 int eN = eN_s[eN_side],
1484 ebN_local = ebN_local_s[eN_side],
1485 kb = kmap[eN_side][k],
1486 eN_nDOF_trial_element = eN * nDOF_trial_element,
1487 ebN_local_kb = ebN_local * nQuadraturePoints_elementBoundary + kb;
1488 double ua = 0.0, ub = 0.0,
1489 grad_ua[nSpace] = {0., 0.}, grad_ub[nSpace] = {0., 0.},
1490 jac_int[nSpace * nSpace], jacDet_int, jacInv_int[nSpace * nSpace],
1491 boundaryJac[nSpace * (nSpace - 1)],
1492 metricTensor[(nSpace - 1) * (nSpace - 1)], metricTensorDetSqrt,
1493 normal[2], x_int, y_int, z_int;
1494 ck.calculateMapping_elementBoundary(eN, ebN_local, kb, ebN_local_kb,
1495 mesh_dof.data(), mesh_l2g.data(),
1496 mesh_trial_trace_ref.data(), mesh_grad_trial_trace_ref.data(),
1497 boundaryJac_ref.data(), jac_int, jacDet_int, jacInv_int,
1498 boundaryJac, metricTensor, metricTensorDetSqrt,
1499 normal_ref.data(), normal, x_int, y_int, z_int);
1500 // dS and the reference normal n_e are fixed from side 0; side 1's own
1501 // outward normal is -n_e, which sign_ne converts back.
1502 if (eN_side == 0)
1503 dS = metricTensorDetSqrt * dS_ref.data()[kb];
1504 double sign_ne = (eN_side == 0) ? 1.0 : -1.0;
1505
1506 auto element_u = xt::pyarray<double>::from_shape({nDOF_trial_element});
1507 double element_phi_f[nDOF_trial_element], element_nodes[nDOF_trial_element * 3];
1508 for (int i = 0; i < nDOF_trial_element; i++)
1509 {
1510 int eN_i = eN * nDOF_trial_element + i;
1511 element_u.data()[i] = u_dof.data()[u_l2g.data()[eN_i]];
1512 element_phi_f[i] = immersedBoundary_sdf_nodes.data()[u_l2g.data()[eN_i]];
1513 for (int I = 0; I < 3; I++)
1514 element_nodes[i * 3 + I] = mesh_dof.data()[u_l2g.data()[eN_i] * 3 + I];
1515 }
1516 if (gf_f_boundary_gen[eN] != ifemGeometryGeneration)
1517 {
1518 gf_f_boundary_icase[eN] = gf_f_cache[eN].calculate(element_phi_f, element_nodes, xB_ref_faces, mua, mub, jf, true, false);
1519 gf_f_boundary_gen[eN] = ifemGeometryGeneration;
1520 }
1521 GfType &gf_f = gf_f_cache[eN];
1522 gf_f.set_boundary_quad(ebN_local_kb);
1523
1524 double va[nDOF_trial_element], va_grad_trial[nDOF_trial_element * nSpace],
1525 vb[nDOF_trial_element], vb_grad_trial[nDOF_trial_element * nSpace];
1526 for (int i = 0; i < nDOF_trial_element; i++)
1527 {
1528 va[i] = gf_f.VA(i);
1529 va_grad_trial[i * nSpace + 0] = gf_f.VA_x(i);
1530 va_grad_trial[i * nSpace + 1] = gf_f.VA_y(i);
1531 vb[i] = gf_f.VB(i);
1532 vb_grad_trial[i * nSpace + 0] = gf_f.VB_x(i);
1533 vb_grad_trial[i * nSpace + 1] = gf_f.VB_y(i);
1534 }
1535 ck.valFromElementDOF(element_u.data(), va, ua);
1536 ck.gradFromElementDOF(element_u.data(), va_grad_trial, grad_ua);
1537 ck.valFromElementDOF(element_u.data(), vb, ub);
1538 ck.gradFromElementDOF(element_u.data(), vb_grad_trial, grad_ub);
1539
1540 ua_s[eN_side] = ua;
1541 ub_s[eN_side] = ub;
1542 double fa = 0.0, fb = 0.0;
1543 for (int I = 0; I < nSpace; I++)
1544 {
1545 fa += mua * grad_ua[I] * normal[I];
1546 fb += mub * grad_ub[I] * normal[I];
1547 }
1548 flux_a_s[eN_side] = sign_ne * fa;
1549 flux_b_s[eN_side] = sign_ne * fb;
1550
1551 for (int i = 0; i < nDOF_test_element; i++)
1552 {
1553 int dof_i = u_l2g.data()[eN_nDOF_trial_element + i];
1554 va_m[eN_side][dof_i] = va[i];
1555 vb_m[eN_side][dof_i] = vb[i];
1556 double fta = 0.0, ftb = 0.0;
1557 for (int I = 0; I < nSpace; I++)
1558 {
1559 fta += mua * va_grad_trial[i * nSpace + I] * normal[I];
1560 ftb += mub * vb_grad_trial[i * nSpace + I] * normal[I];
1561 }
1562 fta_m[eN_side][dof_i] = sign_ne * fta;
1563 ftb_m[eN_side][dof_i] = sign_ne * ftb;
1564 }
1565 } // eN_side
1566
1567 // t of this physical point along side 0's edge direction, then H_hat(t)
1568 double proj = 0.0;
1569 for (int I = 0; I < 3; I++)
1570 proj += (xq[0][k][I] - edge_nodes[0][I]) * edge_vec[I];
1571 double t_edge = (edge_len2 > 1.0e-24) ? (proj / edge_len2) : 0.0;
1572 double H_e = equivalent_polynomials::evaluate_edge_poly<nP_edge>(C_H_edge, t_edge);
1573 double ImH_e = 1.0 - H_e;
1574
1575 double avg_flux_a = 0.5 * (flux_a_s[0] + flux_a_s[1]);
1576 double avg_flux_b = 0.5 * (flux_b_s[0] + flux_b_s[1]);
1577 double jump_ua = ua_s[0] - ua_s[1];
1578 double jump_ub = ub_s[0] - ub_s[1];
1579
1580 // Blend the region-wise PRODUCTS, never the individual factors:
1581 // (ImH*A_f + H*B_f)*(ImH*A_v + H*B_v) != ImH*A_f*A_v + H*B_f*B_v.
1582 for (int side = 0; side < 2; side++)
1583 {
1584 double v_sign = (side == 0) ? 1.0 : -1.0;
1585 for (std::map<int, double>::iterator w = va_m[side].begin(); w != va_m[side].end(); ++w)
1586 {
1587 int dof_i = w->first;
1588 double jva_i = v_sign * va_m[side][dof_i], jvb_i = v_sign * vb_m[side][dof_i];
1589 double Pa = avg_flux_a * jva_i + 0.5 * fta_m[side][dof_i] * jump_ua;
1590 double Pb = avg_flux_b * jvb_i + 0.5 * ftb_m[side][dof_i] * jump_ub;
1591 globalResidual.data()[offset_u + stride_u * dof_i] -=
1592 immersedSCIFEM_switch * (ImH_e * Pa + H_e * Pb) * dS;
1593 // Interior-penalty stabilization gamma/h * int_e [[u]][[v]] ds. eq (3.7) is
1594 // consistency + adjoint only, so it carries no coercivity guarantee; this is
1595 // the standard remedy. Added with a PLUS sign (like the bulk stiffness) so it
1596 // is positive semi-definite, and blended by the same edge Heaviside so it
1597 // stays exact for P1 (integrand degree <= 2) and P2 (degree <= 4 = nP_edge).
1598 globalResidual.data()[offset_u + stride_u * dof_i] +=
1599 (immersedSCIFEM_penalty / h_edge) *
1600 (ImH_e * jump_ua * jva_i + H_e * jump_ub * jvb_i) * dS;
1601 }
1602 }
1603 } // k
1604 ++it;
1605 }
1606 else
1607 {
1608 it = ifem_boundaries.erase(it);
1609 }
1610 } // ifem element boundaries
1611 //
1612 // loop over exterior element boundaries to calculate surface integrals and load into element and global residuals
1613 //
1614 // ebNE is the Exterior element boundary INdex
1615 // ebN is the element boundary INdex
1616 // eN is the element index
1617 for (int ebNE = 0; ebNE < nExteriorElementBoundaries_global; ebNE++)
1618 {
1619 int ebN = exteriorElementBoundariesArray.data()[ebNE],
1620 eN = elementBoundaryElementsArray.data()[ebN * 2 + 0],
1621 ebN_local = elementBoundaryLocalElementBoundariesArray.data()[ebN * 2 + 0],
1622 eN_nDOF_trial_element = eN * nDOF_trial_element;
1623 double elementResidual_u[nDOF_test_element];
1624 for (int i = 0; i < nDOF_test_element; i++)
1625 {
1626 elementResidual_u[i] = 0.0;
1627 }
1628 for (int kb = 0; kb < nQuadraturePoints_elementBoundary; kb++)
1629 {
1630 int ebNE_kb = ebNE * nQuadraturePoints_elementBoundary + kb,
1631 ebNE_kb_nSpace = ebNE_kb * nSpace,
1632 ebN_local_kb = ebN_local * nQuadraturePoints_elementBoundary + kb,
1633 ebN_local_kb_nSpace = ebN_local_kb * nSpace;
1634 double u_ext = 0.0,
1635 grad_u_ext[nSpace],
1636 m_ext = 0.0,
1637 dm_ext = 0.0,
1638 *a_ext,
1639 /* *da_exxt, */
1640 f_ext[nSpace],
1641 df_ext[nSpace],
1642 r_ext = 0.0,
1643 /* dr_ext=0.0, */
1644 flux_diff_ext = 0.0,
1645 flux_advect_ext = 0.0,
1646 bc_u_ext = 0.0,
1647 // bc_grad_u_ext[nSpace],
1648 bc_m_ext = 0.0,
1649 bc_dm_ext = 0.0,
1650 bc_f_ext[nSpace],
1651 bc_df_ext[nSpace],
1652 jac_ext[nSpace * nSpace],
1653 jacDet_ext,
1654 jacInv_ext[nSpace * nSpace],
1655 boundaryJac[nSpace * (nSpace - 1)],
1656 metricTensor[(nSpace - 1) * (nSpace - 1)],
1657 metricTensorDetSqrt,
1658 dS,
1659 u_test_dS[nDOF_test_element],
1660 u_grad_trial_trace[nDOF_trial_element * nSpace],
1661 u_grad_test_dS[nDOF_trial_element * nSpace],
1662 normal[nSpace], x_ext, y_ext, z_ext, xt_ext, yt_ext, zt_ext, integralScaling,
1663 //
1664 G[nSpace * nSpace], G_dd_G, tr_G;
1665 //
1666 // calculate the solution and gradients at quadrature points
1667 //
1668 // compute information about mapping from reference element to physical element
1669 ck.calculateMapping_elementBoundary(eN,
1670 ebN_local,
1671 kb,
1672 ebN_local_kb,
1673 mesh_dof.data(),
1674 mesh_l2g.data(),
1675 mesh_trial_trace_ref.data(),
1676 mesh_grad_trial_trace_ref.data(),
1677 boundaryJac_ref.data(),
1678 jac_ext,
1679 jacDet_ext,
1680 jacInv_ext,
1681 boundaryJac,
1682 metricTensor,
1683 metricTensorDetSqrt,
1684 normal_ref.data(),
1685 normal,
1686 x_ext, y_ext, z_ext);
1687 /* ck.calculateMappingVelocity_elementBoundary(eN, */
1688 /* ebN_local, */
1689 /* kb, */
1690 /* ebN_local_kb, */
1691 /* mesh_velocity_dof, */
1692 /* mesh_l2g, */
1693 /* mesh_trial_trace_ref, */
1694 /* xt_ext,yt_ext,zt_ext, */
1695 /* normal, */
1696 /* boundaryJac, */
1697 /* metricTensor, */
1698 /* integralScaling); */
1699 /*dS = ((1.0-MOVING_DOMAIN)*metricTensorDetSqrt + MOVING_DOMAIN*integralScaling)*dS_ref.data()[kb];*/
1700 dS = metricTensorDetSqrt * dS_ref.data()[kb];
1701 // get the metric tensor
1702 // cek todo use symmetry
1703 ck.calculateG(jacInv_ext, G, G_dd_G, tr_G);
1704 // compute shape and solution information
1705 // shape
1706 // std::cout << "Calculating gradTrialFromRef from calculateResidual() 3" << std::endl;
1707 ck.gradTrialFromRef(&u_grad_trial_trace_ref.data()[ebN_local_kb_nSpace * nDOF_trial_element], jacInv_ext, u_grad_trial_trace);
1708 // solution and gradients
1709 ck.valFromDOF(u_dof.data(), &u_l2g.data()[eN_nDOF_trial_element], &u_trial_trace_ref.data()[ebN_local_kb * nDOF_test_element], u_ext);
1710 ck.gradFromDOF(u_dof.data(), &u_l2g.data()[eN_nDOF_trial_element], u_grad_trial_trace, grad_u_ext);
1711 // precalculate test function products with integration weights
1712 for (int j = 0; j < nDOF_trial_element; j++)
1713 {
1714 u_test_dS[j] = u_test_trace_ref.data()[ebN_local_kb * nDOF_test_element + j] * dS;
1715 for (int I = 0; I < nSpace; I++)
1716 u_grad_test_dS[j * nSpace + I] = u_grad_trial_trace[j * nSpace + I] * dS; // cek hack, using trial
1717 }
1718 //
1719 // load the boundary values
1720 //
1721 bc_u_ext = isDOFBoundary_u.data()[ebNE_kb] * ebqe_bc_u_ext.data()[ebNE_kb] + (1 - isDOFBoundary_u.data()[ebNE_kb]) * u_ext;
1722 //
1723 //
1724 // calculate the pde coefficients using the solution and the boundary values for the solution
1725 //
1726 a_ext = &ebqe_a.data()[ebNE_kb * sd_rowptr[nSpace]];
1727 for (int I = 0; I < nSpace; I++)
1728 {
1729 f_ext[I] = ebqe_v.data()[ebNE_kb * nSpace + I] * u_ext;
1730 df_ext[I] = ebqe_v.data()[ebNE_kb * nSpace + I];
1731 bc_f_ext[I] = ebqe_v.data()[ebNE_kb * nSpace + I] * bc_u_ext;
1732 bc_df_ext[I] = ebqe_v.data()[ebNE_kb * nSpace + I];
1733 }
1734 /* evaluateCoefficients(&ebqe_velocity_ext.data()[ebNE_kb_nSpace], */
1735 /* u_ext, */
1736 /* //VRANS */
1737 /* porosity_ext, */
1738 /* // */
1739 /* m_ext, */
1740 /* dm_ext, */
1741 /* f_ext, */
1742 /* df_ext); */
1743 /* evaluateCoefficients(&ebqe_velocity_ext.data()[ebNE_kb_nSpace], */
1744 /* bc_u_ext, */
1745 /* //VRANS */
1746 /* porosity_ext, */
1747 /* // */
1748 /* bc_m_ext, */
1749 /* bc_dm_ext, */
1750 /* bc_f_ext, */
1751 /* bc_df_ext); */
1752 //
1753 // moving mesh
1754 //
1755 /* double velocity_ext[nSpace]; */
1756 /* double mesh_velocity[3]; */
1757 /* mesh_velocity[0] = xt_ext; */
1758 /* mesh_velocity[1] = yt_ext; */
1759 /* mesh_velocity[2] = zt_ext; */
1760 /* for (int I=0;I<nSpace;I++) */
1761 /* velocity_ext[I] = ebqe_velocity_ext.data()[ebNE_kb_nSpace+0] - MOVING_DOMAIN*mesh_velocity[I]; */
1762 //
1763 // calculate the numerical fluxes
1764 //
1765 exteriorNumericalDiffusiveFlux(sd_rowptr.data(),
1766 sd_colind.data(),
1767 isDOFBoundary_u.data()[ebNE_kb],
1768 isDiffusiveFluxBoundary_u.data()[ebNE_kb],
1769 normal,
1770 a_ext,
1771 bc_u_ext,
1772 ebqe_bc_flux_u_ext.data()[ebNE_kb],
1773 a_ext,
1774 grad_u_ext,
1775 u_ext,
1776 ebqe_penalty_ext.data()[ebNE_kb],
1777 flux_diff_ext);
1778 exteriorNumericalAdvectiveFlux(isDOFBoundary_u.data()[ebNE_kb],
1779 isAdvectiveFluxBoundary_u.data()[ebNE_kb],
1780 normal,
1781 bc_u_ext,
1782 ebqe_bc_flux_u_ext.data()[ebNE_kb],
1783 u_ext,
1784 df_ext,
1785 flux_advect_ext);
1786 // ebqe_flux.data()[ebNE_kb] = flux_ext;
1787 //
1788 // update residuals
1789 //
1790 for (int i = 0; i < nDOF_test_element; i++)
1791 {
1792 // int ebNE_kb_i = ebNE_kb*nDOF_test_element+i;
1793 elementResidual_u[i] += ck.ExteriorElementBoundaryFlux(flux_diff_ext + flux_advect_ext, u_test_dS[i]) +
1794 ck.ExteriorElementBoundaryDiffusionAdjoint(isDOFBoundary_u.data()[ebNE_kb],
1795 isDiffusiveFluxBoundary_u.data()[ebNE_kb],
1796 eb_adjoint_sigma,
1797 u_ext,
1798 bc_u_ext,
1799 normal,
1800 sd_rowptr.data(),
1801 sd_colind.data(),
1802 a_ext,
1803 &u_grad_test_dS[i * nSpace]);
1804 } // i
1805 } // kb
1806 //
1807 // update the element and global residual storage
1808 //
1809 for (int i = 0; i < nDOF_test_element; i++)
1810 {
1811 int eN_i = eN * nDOF_test_element + i;
1812 globalResidual.data()[offset_u + stride_u * u_l2g.data()[eN_i]] += elementResidual_u[i];
1813 } // i
1814 } // ebNE
1815 }
1816
1817 inline void calculateElementJacobian(int icase_f,
1818 // element
1819 xt::pyarray<double> &mesh_trial_ref,
1820 xt::pyarray<double> &mesh_grad_trial_ref,
1821 xt::pyarray<double> &mesh_dof,
1822 xt::pyarray<int> &mesh_l2g,
1823 xt::pyarray<double> &x_ref,
1824 xt::pyarray<double> &dV_ref,
1825 xt::pyarray<double> &u_trial_ref,
1826 xt::pyarray<double> &u_grad_trial_ref,
1827 xt::pyarray<double> &u_test_ref,
1828 xt::pyarray<double> &u_grad_test_ref,
1829 xt::pyarray<double> &elementDiameter,
1830 xt::pyarray<double> &elementBoundaryDiameter,
1831 xt::pyarray<double> &nodeDiametersArray,
1832 xt::pyarray<double> &cfl,
1833 double Ct_sge,
1834 double sc_uref,
1835 double sc_alpha,
1836 double useMetrics,
1837 // element boundary
1838 xt::pyarray<double> &mesh_trial_trace_ref,
1839 xt::pyarray<double> &mesh_grad_trial_trace_ref,
1840 xt::pyarray<double> &dS_ref,
1841 xt::pyarray<double> &u_trial_trace_ref,
1842 xt::pyarray<double> &u_grad_trial_trace_ref,
1843 xt::pyarray<double> &u_test_trace_ref,
1844 xt::pyarray<double> &u_grad_test_trace_ref,
1845 xt::pyarray<double> &normal_ref,
1846 xt::pyarray<double> &boundaryJac_ref,
1847 // physics
1848 int nElements_global,
1849 int nElementBoundaries_owned,
1850 xt::pyarray<int> &u_l2g,
1851 xt::pyarray<double> &u_dof,
1852 xt::pyarray<int> &sd_rowptr,
1853 xt::pyarray<int> &sd_colind,
1854 xt::pyarray<double> &q_a,
1855 xt::pyarray<double> &q_v,
1856 xt::pyarray<double> &q_r,
1857 int lag_shockCapturing,
1858 double shockCapturingDiffusion,
1859 xt::pyarray<double> &q_numDiff_u,
1860 xt::pyarray<double> &q_numDiff_u_last,
1861 xt::pyarray<double> &elementJacobian_u_u,
1862 xt::pyarray<double> &element_u,
1863 int eN,
1864 const bool embeddedBoundary,
1865 const double embeddedBoundary_penalty,
1866 xt::pyarray<double> &embeddedBoundary_normal_q,
1867 xt::pyarray<double> &embeddedBoundary_u_q,
1868 const bool immersedBoundary,
1869 const double immersedBoundary_penalty,
1870 xt::pyarray<double> &immersedBoundary_sdf_q,
1871 xt::pyarray<double> &immersedBoundary_normal_q,
1872 xt::pyarray<double> &immersedBoundary_u_q,
1873 xt::pyarray<double> &immersedBoundary_fluxJump_q,
1874 xt::pyarray<double> &immersedBoundary_fluxJumpVector_q,
1875 xt::pyarray<double> &immersedBoundary_solutionJump_nodes,
1876 double *element_phi_f,
1877 double mua,
1878 double mub,
1879 const bool PG)
1880 {
1881 // per-element cached equivalent-polynomial/IFEM reconstruction
1882 // (see ensureIFEMCacheSized / ifemGeometryGeneration)
1883 GfType &gf_f = gf_f_cache[eN];
1884 GfType &gf_s = gf_s_cache[eN];
1885 // std::cout << "Calculating element Jacobian for element " << eN << std::endl;
1886 for (int i = 0; i < nDOF_test_element; i++)
1887 for (int j = 0; j < nDOF_trial_element; j++)
1888 {
1889 elementJacobian_u_u.data()[i * nDOF_trial_element + j] = 0.0;
1890 }
1891 for (int k = 0; k < nQuadraturePoints_element; k++)
1892 {
1893 // std::cout << " quadrature point " << k << std::endl;
1894 gf_s.set_quad(k);
1895 gf_f.set_quad(k);
1896 int eN_k = eN * nQuadraturePoints_element + k; // index to a scalar at a quadrature point
1897 int eN_k_3d = eN_k * 3;
1898 // declare local storage
1899 double u = 0.0,
1900 grad_u[nSpace],
1901 ua = 0.0,
1902 grad_ua[nSpace],
1903 ub = 0.0,
1904 grad_ub[nSpace],
1905 m = 0.0, dm = 0.0,
1906 h_phi = 0.0,
1907 r_s = 0.0, dr_s = 0.0,
1908 f[nSpace], df[nSpace],
1909 f_s[nSpace] = {0., 0.}, df_s[nSpace] = {0., 0.},
1910 ham_s = 0.0, dham_s[nSpace] = {0., 0.},
1911 r_f = 0.0, dr_f = 0.0,
1912 f_f[nSpace] = {0., 0.}, df_f[nSpace] = {0., 0.},
1913 ham_f = 0.0, dham_f[nSpace] = {0., 0.},
1914 m_t = 0.0, dm_t = 0.0,
1915 dpdeResidual_u_u[nDOF_trial_element],
1916 Lstar_u[nDOF_test_element],
1917 dsubgridError_u_u[nDOF_trial_element],
1918 tau = 0.0, tau0 = 0.0, tau1 = 0.0,
1919 *a = NULL,
1920 dr = 0.0,
1921 jac[nSpace * nSpace],
1922 jacDet,
1923 jacInv[nSpace * nSpace],
1924 u_grad_trial[nDOF_trial_element * nSpace],
1925 ua_grad_trial[nDOF_trial_element * nSpace],
1926 ub_grad_trial[nDOF_trial_element * nSpace],
1927 dV,
1928 ua_trial[nDOF_trial_element],
1929 ub_trial[nDOF_trial_element],
1930 u_test_dV[nDOF_test_element],
1931 ua_test_dV[nDOF_test_element],
1932 ub_test_dV[nDOF_test_element],
1933 u_grad_test_dV[nDOF_test_element * nSpace],
1934 ua_grad_test_dV[nDOF_test_element * nSpace],
1935 ub_grad_test_dV[nDOF_test_element * nSpace],
1936 x, y, z,
1937 G[nSpace * nSpace], G_dd_G, tr_G;
1938 //
1939 // calculate solution and gradients at quadrature points
1940 //
1941 ck.calculateMapping_element(eN,
1942 k,
1943 mesh_dof.data(),
1944 mesh_l2g.data(),
1945 mesh_trial_ref.data(),
1946 mesh_grad_trial_ref.data(),
1947 jac,
1948 jacDet,
1949 jacInv,
1950 x, y, z);
1951 ck.calculateH_element(eN,
1952 k,
1953 nodeDiametersArray.data(),
1954 mesh_l2g.data(),
1955 mesh_trial_ref.data(),
1956 h_phi);
1957 // get the physical integration weight
1958 dV = fabs(jacDet) * dV_ref.data()[k];
1959 // get metric tensor and friends
1960 ck.calculateG(jacInv, G, G_dd_G, tr_G);
1961 // get the trial function gradients
1962 // std::cout << " calculating gradTrialfromRef from calculateElementJacobian() "<< std::endl;
1963 ck.gradTrialFromRef(&u_grad_trial_ref.data()[k * nDOF_trial_element * nSpace], jacInv, u_grad_trial);
1964 // get the solution
1965 ck.valFromElementDOF(element_u.data(), &u_trial_ref.data()[k * nDOF_trial_element], u);
1966 // get the solution gradients
1967 ck.gradFromElementDOF(element_u.data(), u_grad_trial, grad_u);
1968 // precalculate test function products with integration weights
1969 for (int j = 0; j < nDOF_trial_element; j++)
1970 {
1971 u_test_dV[j] = u_test_ref.data()[k * nDOF_trial_element + j] * dV;
1972 for (int I = 0; I < nSpace; I++)
1973 {
1974 u_grad_test_dV[j * nSpace + I] = u_grad_trial[j * nSpace + I] * dV; // cek warning won't work for Petrov-Galerkin
1975 //PGIFEM
1976 //ua_grad_test_dV[j * nSpace + I] = u_grad_trial[j * nSpace + I] * dV; // cek warning won't work for Petrov-Galerkin
1977 //ub_grad_test_dV[j * nSpace + I] = u_grad_trial[j * nSpace + I] * dV; // cek warning won't work for Petrov-Galerkin
1978 }
1979 }
1980 if (icase_f == 0)
1981 {
1982 double va[nDOF_trial_element], va_grad_trial[nDOF_trial_element * nSpace],
1983 vb[nDOF_trial_element], vb_grad_trial[nDOF_trial_element * nSpace];
1984 for (int i = 0; i < nDOF_trial_element; i++)
1985 {
1986 va[i] = gf_f.VA(i);
1987 va_grad_trial[i * nSpace + 0] = gf_f.VA_x(i);
1988 va_grad_trial[i * nSpace + 1] = gf_f.VA_y(i);
1989 vb[i] = gf_f.VB(i);
1990 vb_grad_trial[i * nSpace + 0] = gf_f.VB_x(i);
1991 vb_grad_trial[i * nSpace + 1] = gf_f.VB_y(i);
1992 // std::cout << " va[" << i << "]=" << va[i] << std::endl;
1993 // std::cout << " vb[" << i << "]=" << vb[i] << std::endl;
1994 // std::cout << " va_grad_trial[" << i << "]=" << va_grad_trial[i * nSpace + 0] << "," << va_grad_trial[i * nSpace + 1] << std::endl;
1995 // std::cout << " vb_grad_trial[" << i << "]=" << vb_grad_trial[i * nSpace + 0] << "," << vb_grad_trial[i * nSpace + 1] << std::endl;
1996 }
1997 ck.valFromElementDOF(element_u.data(), va, ua);
1998 ck.gradFromElementDOF(element_u.data(), va_grad_trial, grad_ua);
1999 ck.valFromElementDOF(element_u.data(), vb, ub);
2000 ck.gradFromElementDOF(element_u.data(), vb_grad_trial, grad_ub);
2001 for (int i = 0; i < nDOF_test_element; i++)
2002 {
2003 ua_test_dV[i] = va[i] * dV;
2004 ub_test_dV[i] = vb[i] * dV;
2005 for (int I = 0; I < nSpace; I++)
2006 {
2007 ua_grad_test_dV[i * nSpace + I] = va_grad_trial[i * nSpace + I] * dV;
2008 ub_grad_test_dV[i * nSpace + I] = vb_grad_trial[i * nSpace + I] * dV;
2009 }
2010 }
2011 for (int j = 0; j < nDOF_trial_element; j++)
2012 {
2013 ua_trial[j] = va[j];
2014 ub_trial[j] = vb[j];
2015 for (int I = 0; I < nSpace; I++)
2016 {
2017 ua_grad_trial[j * nSpace + I] = va_grad_trial[j * nSpace + I];
2018 ub_grad_trial[j * nSpace + I] = vb_grad_trial[j * nSpace + I];
2019 }
2020 // std::cout << " ua_grad_test_dV[" << j << "]=" << ua_grad_test_dV[j * nSpace + 0] << "," << ua_grad_test_dV[j * nSpace + 1] << std::endl;
2021 // std::cout << " ub_grad_test_dV[" << j << "]=" << ub_grad_test_dV[j * nSpace + 0] << "," << ub_grad_test_dV[j * nSpace + 1] << std::endl;
2022 }
2023 if (PG)
2024 {
2025 // Same substitution as the residual; trial-side (j-index) arrays
2026 // are untouched, only the test side (i-index) changes.
2027 for (int i = 0; i < nDOF_test_element; i++)
2028 {
2029 ua_test_dV[i] = u_test_dV[i];
2030 ub_test_dV[i] = u_test_dV[i];
2031 for (int I = 0; I < nSpace; I++)
2032 {
2033 ua_grad_test_dV[i * nSpace + I] = u_grad_test_dV[i * nSpace + I];
2034 ub_grad_test_dV[i * nSpace + I] = u_grad_test_dV[i * nSpace + I];
2035 }
2036 }
2037 }
2038 }
2039 //
2040 // calculate pde coefficients and derivatives at quadrature points
2041 //
2042 // evaluateCoefficients()
2043 a = &q_a.data()[eN_k * sd_rowptr.data()[nSpace]];
2044 for (int I = 0; I < nSpace; I++)
2045 df[I] = q_v.data()[eN_k * nSpace + I];
2046 dr = 0.0;
2047 const double H_s = gf_s.H(0., 0.);
2048 const double D_s = gf_s.D(0., 0.);
2049 if (embeddedBoundary)
2050 {
2051 double level_set_normal[nSpace];
2052 double sign = 0.0;
2053 double norm_exact = 0.0, norm_cut = 0.0;
2054 for (int I = 0; I < nSpace; I++)
2055 {
2056 sign += embeddedBoundary_normal_q.data()[eN_k_3d + I] * gf_s.get_normal()[I];
2057 level_set_normal[I] = gf_s.get_normal()[I];
2058 norm_cut += level_set_normal[I] * level_set_normal[I];
2059 norm_exact += embeddedBoundary_normal_q.data()[eN_k_3d + I] * embeddedBoundary_normal_q.data()[eN_k_3d + I];
2060 }
2061 assert(std::fabs(1.0 - norm_cut) < 1.0e-8);
2062 assert(std::fabs(1.0 - norm_exact) < 1.0e-8);
2063 if (sign < 0.0)
2064 for (int I = 0; I < nSpace; I++)
2065 level_set_normal[I] *= -1.0;
2066 updateEmbeddedBoundaryTerms(embeddedBoundary_penalty / h_phi, // penalty,
2067 dV,
2068 level_set_normal,
2069 embeddedBoundary_u_q.data()[eN_k],
2070 u,
2071 grad_u,
2072 a[0], // assume scalar diffusion for now
2073 r_s,
2074 dr_s,
2075 ham_s,
2076 dham_s,
2077 f_s,
2078 df_s,
2079 D_s);
2080 }
2081 const double ImH_f = gf_f.ImH(0., 0.);
2082 const double H_f = gf_f.H(0., 0.);
2083 const double D_f = gf_f.D(0., 0.);
2084 if (immersedBoundary)
2085 {
2086 double level_set_normal[nSpace];
2087 double sign = 0.0;
2088 double norm_exact = 0.0, norm_cut = 0.0;
2089 for (int I = 0; I < nSpace; I++)
2090 {
2091 sign += immersedBoundary_normal_q.data()[eN_k_3d + I] * gf_f.get_normal()[I];
2092 level_set_normal[I] = gf_f.get_normal()[I];
2093 norm_cut += level_set_normal[I] * level_set_normal[I];
2094 norm_exact += immersedBoundary_normal_q.data()[eN_k_3d + I] * immersedBoundary_normal_q.data()[eN_k_3d + I];
2095 }
2096 assert(std::fabs(1.0 - norm_cut) < 1.0e-8);
2097 assert(std::fabs(1.0 - norm_exact) < 1.0e-8);
2098 if (sign < 0.0)
2099 for (int I = 0; I < nSpace; I++)
2100 level_set_normal[I] *= -1.0;
2101 updateImmersedBoundaryTerms(immersedBoundary_penalty / h_phi, // penalty,
2102 dV,
2103 level_set_normal,
2104 x,
2105 y,
2106 z,
2107 immersedBoundary_u_q.data()[eN_k],
2108 u,
2109 grad_u,
2110 a[0], // assume scalar diffusion for now
2111 r_f,
2112 dr_f,
2113 ham_f,
2114 dham_f,
2115 f_f,
2116 df_f,
2117 immersedBoundary_fluxJump_q.data()[eN_k],
2118 &immersedBoundary_fluxJumpVector_q.data()[eN_k_3d],
2119 D_f);
2120 }
2121 //
2122 // calculate subgrid error contribution to the Jacobian (strong residual, adjoint, jacobian of strong residual)
2123 //
2124 // calculate the adjoint times the test functions
2125 for (int i = 0; i < nDOF_test_element; i++)
2126 {
2127 // int eN_k_i_nSpace = (eN_k*nDOF_trial_element+i)*nSpace;
2128 // Lstar_u[i]=ck.Advection_adjoint(df,&u_grad_test_dV.data()[eN_k_i_nSpace]);
2129 int i_nSpace = i * nSpace;
2130 Lstar_u[i] = ck.Advection_adjoint(df, &u_grad_test_dV[i_nSpace]);
2131 }
2132 // calculate the Jacobian of strong residual
2133 for (int j = 0; j < nDOF_trial_element; j++)
2134 {
2135 // int eN_k_j=eN_k*nDOF_trial_element+j;
2136 // int eN_k_j_nSpace = eN_k_j*nSpace;
2137 int j_nSpace = j * nSpace;
2138 dpdeResidual_u_u[j] = ck.MassJacobian_strong(dm_t, u_trial_ref.data()[k * nDOF_trial_element + j]) +
2139 ck.AdvectionJacobian_strong(df, &u_grad_trial[j_nSpace]);
2140 }
2141 // tau and tau*Res
2142 calculateSubgridError_tau(elementDiameter.data()[eN],
2143 dm_t,
2144 df,
2145 cfl.data()[eN_k],
2146 tau0);
2147
2148 calculateSubgridError_tau(Ct_sge,
2149 G,
2150 dm_t,
2151 df,
2152 tau1,
2153 cfl.data()[eN_k]);
2154 tau = useMetrics * tau1 + (1.0 - useMetrics) * tau0;
2155
2156 for (int j = 0; j < nDOF_trial_element; j++)
2157 dsubgridError_u_u[j] = -tau * dpdeResidual_u_u[j];
2158
2159 // Leveque & Li 1994, Examples 1, 3, 4, PWC, PWL, PWQ
2160 double a_loc[nSpace * nSpace];
2161 for (int I = 0; I < nSpace * nSpace; I++) a_loc[I] = 0.0;
2162 for (int i = 0; i < nDOF_test_element; i++)
2163 {
2164 int i_nSpace = i * nSpace;
2165 for (int j = 0; j < nDOF_trial_element; j++)
2166 {
2167 int j_nSpace = j * nSpace;
2168 if (icase_f == 0)
2169 {
2170 if(!gf_f.exact.edge && !gf_f.exact.corner)
2171 {
2172 for (int I = 0; I < nSpace; I++) a_loc[I * nSpace + I] = mua;
2173 elementJacobian_u_u.data()[i * nDOF_trial_element + j] += ImH_f * H_s * (ck.AdvectionJacobian_weak(df, ua_trial[j], &ua_grad_test_dV[i_nSpace]) +
2174 ck.SimpleDiffusionJacobian_weak(sd_rowptr.data(), sd_colind.data(), a_loc, &ua_grad_trial[j_nSpace], &ua_grad_test_dV[i_nSpace]) +
2175 ck.ReactionJacobian_weak(dr, ua_trial[j], ua_test_dV[i]) +
2176 ck.NumericalDiffusionJacobian(q_numDiff_u_last.data()[eN_k], &ua_grad_trial[j_nSpace], &ua_grad_test_dV[i_nSpace]));
2177
2178 for (int I = 0; I < nSpace; I++) a_loc[I * nSpace + I] = mub;
2179 elementJacobian_u_u.data()[i * nDOF_trial_element + j] += H_f * H_s * (ck.AdvectionJacobian_weak(df, ub_trial[j], &ub_grad_test_dV[i_nSpace]) +
2180 ck.SimpleDiffusionJacobian_weak(sd_rowptr.data(), sd_colind.data(), a_loc, &ub_grad_trial[j_nSpace], &ub_grad_test_dV[i_nSpace]) +
2181 ck.ReactionJacobian_weak(dr, ub_trial[j], ub_test_dV[i]) +
2182 ck.NumericalDiffusionJacobian(q_numDiff_u_last.data()[eN_k], &ub_grad_trial[j_nSpace], &ub_grad_test_dV[i_nSpace]));
2183 }
2184 else if (gf_f.exact.edge == -1 || gf_f.exact.corner == -1)
2185 {
2186 for (int I = 0; I < nSpace; I++) a_loc[I * nSpace + I] = mua;
2187 elementJacobian_u_u.data()[i * nDOF_trial_element + j] += ImH_f * H_s * (ck.AdvectionJacobian_weak(df, ua_trial[j], &ua_grad_test_dV[i_nSpace]) +
2188 ck.SimpleDiffusionJacobian_weak(sd_rowptr.data(), sd_colind.data(), a_loc, &ua_grad_trial[j_nSpace], &ua_grad_test_dV[i_nSpace]) +
2189 ck.ReactionJacobian_weak(dr, ua_trial[j], ua_test_dV[i]) +
2190 ck.NumericalDiffusionJacobian(q_numDiff_u_last.data()[eN_k], &ua_grad_trial[j_nSpace], &ua_grad_test_dV[i_nSpace]));
2191 }
2192 else if (gf_f.exact.edge == 1 || gf_f.exact.corner == 1)
2193 {
2194 for (int I = 0; I < nSpace; I++) a_loc[I * nSpace + I] = mub;
2195 elementJacobian_u_u.data()[i * nDOF_trial_element + j] += H_f * H_s * (ck.AdvectionJacobian_weak(df, ub_trial[j], &ub_grad_test_dV[i_nSpace]) +
2196 ck.SimpleDiffusionJacobian_weak(sd_rowptr.data(), sd_colind.data(), a_loc, &ub_grad_trial[j_nSpace], &ub_grad_test_dV[i_nSpace]) +
2197 ck.ReactionJacobian_weak(dr, ub_trial[j], ub_test_dV[i]) +
2198 ck.NumericalDiffusionJacobian(q_numDiff_u_last.data()[eN_k], &ub_grad_trial[j_nSpace], &ub_grad_test_dV[i_nSpace]));
2199 }
2200 else assert(false && "Invalid gf_f.exact.edge/corner values. Should be -1, 0 or +1.");
2201 }
2202 else
2203 {
2204 elementJacobian_u_u.data()[i * nDOF_trial_element + j] += H_s * (ck.AdvectionJacobian_weak(df, u_trial_ref.data()[k * nDOF_trial_element + j], &u_grad_test_dV[i_nSpace]) +
2205 ck.SimpleDiffusionJacobian_weak(sd_rowptr.data(), sd_colind.data(), a, &u_grad_trial[j_nSpace], &u_grad_test_dV[i_nSpace]) +
2206 ck.ReactionJacobian_weak(dr, u_trial_ref.data()[k * nDOF_trial_element + j], u_test_dV[i]) +
2207 ck.SubgridErrorJacobian(dsubgridError_u_u[j], Lstar_u[i]) +
2208 ck.NumericalDiffusionJacobian(q_numDiff_u_last.data()[eN_k], &u_grad_trial[j_nSpace], &u_grad_test_dV[i_nSpace]));
2209 }
2210 if (embeddedBoundary)
2211 {
2212 if (gf_s.exact.edge >=0 && !gf_s.exact.corner)
2213 {
2214 elementJacobian_u_u.data()[i * nDOF_trial_element + j] += (ck.AdvectionJacobian_weak(df_s, u_trial_ref.data()[k * nDOF_trial_element + j], &u_grad_test_dV[i_nSpace])
2215 + ck.ReactionJacobian_weak(dr_s, u_trial_ref.data()[k * nDOF_trial_element + j], u_test_dV[i])
2216 + ck.HamiltonianJacobian_weak(dham_s, &u_grad_trial[j_nSpace], u_test_dV[i]));
2217 }
2218 }
2219 if (immersedBoundary)
2220 {
2221 if (gf_f.exact.edge >=0 && !gf_f.exact.corner){
2222 // std::cout << " last loop: elementJacobian_u_u[" << i << "," << j << "] = " << elementJacobian_u_u.data()[i * nDOF_trial_element + j];
2223 elementJacobian_u_u.data()[i * nDOF_trial_element + j] += (ck.AdvectionJacobian_weak(df_f, u_trial_ref.data()[k * nDOF_trial_element + j], &u_grad_test_dV[i_nSpace]) +
2224 ck.ReactionJacobian_weak(dr_f, u_trial_ref.data()[k * nDOF_trial_element + j], u_test_dV[i]) +
2225 ck.HamiltonianJacobian_weak(dham_f, &u_grad_trial[j_nSpace], u_test_dV[i]));
2226 }
2227 }
2228 // std::cout << " elementJacobian_u_u[" << i << "," << j << "] = " << elementJacobian_u_u.data()[i * nDOF_trial_element + j];
2229 } // j
2230 // std::cout << std::endl;
2231 } // i
2232 // std::cout << std::endl;
2233 } // k
2234 }
2235
2236 void calculateJacobian(arguments_dict &args)
2237 {
2238 xt::pyarray<double> &mesh_trial_ref = args.array<double>("mesh_trial_ref");
2239 xt::pyarray<double> &mesh_grad_trial_ref = args.array<double>("mesh_grad_trial_ref");
2240 xt::pyarray<double> &mesh_dof = args.array<double>("mesh_dof");
2241 xt::pyarray<int> &mesh_l2g = args.array<int>("mesh_l2g");
2242 xt::pyarray<double> &dV_ref = args.array<double>("dV_ref");
2243 xt::pyarray<double> &u_trial_ref = args.array<double>("u_trial_ref");
2244 xt::pyarray<double> &u_grad_trial_ref = args.array<double>("u_grad_trial_ref");
2245 xt::pyarray<double> &u_test_ref = args.array<double>("u_test_ref");
2246 xt::pyarray<double> &u_grad_test_ref = args.array<double>("u_grad_test_ref");
2247 xt::pyarray<double> &elementDiameter = args.array<double>("elementDiameter");
2248 xt::pyarray<double> &cfl = args.array<double>("cfl");
2249 double Ct_sge = args.scalar<double>("Ct_sge");
2250 double sc_uref = args.scalar<double>("sc_uref");
2251 double sc_alpha = args.scalar<double>("sc_alpha");
2252 double useMetrics = args.scalar<double>("useMetrics");
2253 xt::pyarray<double> &mesh_trial_trace_ref = args.array<double>("mesh_trial_trace_ref");
2254 xt::pyarray<double> &mesh_grad_trial_trace_ref = args.array<double>("mesh_grad_trial_trace_ref");
2255 xt::pyarray<double> &dS_ref = args.array<double>("dS_ref");
2256 xt::pyarray<double> &u_trial_trace_ref = args.array<double>("u_trial_trace_ref");
2257 xt::pyarray<double> &u_grad_trial_trace_ref = args.array<double>("u_grad_trial_trace_ref");
2258 xt::pyarray<double> &u_test_trace_ref = args.array<double>("u_test_trace_ref");
2259 xt::pyarray<double> &u_grad_test_trace_ref = args.array<double>("u_grad_test_trace_ref");
2260 xt::pyarray<double> &normal_ref = args.array<double>("normal_ref");
2261 xt::pyarray<double> &boundaryJac_ref = args.array<double>("boundaryJac_ref");
2262 int nElements_global = args.scalar<int>("nElements_global");
2263 xt::pyarray<int> &u_l2g = args.array<int>("u_l2g");
2264 xt::pyarray<double> &u_dof = args.array<double>("u_dof");
2265 xt::pyarray<int> &sd_rowptr = args.array<int>("sd_rowptr");
2266 xt::pyarray<int> &sd_colind = args.array<int>("sd_colind");
2267 xt::pyarray<double> &q_a = args.array<double>("q_a");
2268 xt::pyarray<double> &q_v = args.array<double>("q_v");
2269 xt::pyarray<double> &q_r = args.array<double>("q_r");
2270 int lag_shockCapturing = args.scalar<int>("lag_shockCapturing");
2271 double shockCapturingDiffusion = args.scalar<double>("shockCapturingDiffusion");
2272 xt::pyarray<double> &q_numDiff_u = args.array<double>("q_numDiff_u");
2273 xt::pyarray<double> &q_numDiff_u_last = args.array<double>("q_numDiff_u_last");
2274 xt::pyarray<int> &csrRowIndeces_u_u = args.array<int>("csrRowIndeces_u_u");
2275 xt::pyarray<int> &csrColumnOffsets_u_u = args.array<int>("csrColumnOffsets_u_u");
2276 xt::pyarray<double> &globalJacobian = args.array<double>("globalJacobian");
2277 int nExteriorElementBoundaries_global = args.scalar<int>("nExteriorElementBoundaries_global");
2278 xt::pyarray<int> &exteriorElementBoundariesArray = args.array<int>("exteriorElementBoundariesArray");
2279 xt::pyarray<int> &elementBoundaryElementsArray = args.array<int>("elementBoundaryElementsArray");
2280 xt::pyarray<int> &elementBoundaryLocalElementBoundariesArray = args.array<int>("elementBoundaryLocalElementBoundariesArray");
2281 xt::pyarray<double> &ebqe_a = args.array<double>("ebqe_a");
2282 xt::pyarray<double> &ebqe_v = args.array<double>("ebqe_v");
2283 xt::pyarray<int> &isDOFBoundary_u = args.array<int>("isDOFBoundary_u");
2284 xt::pyarray<double> &ebqe_bc_u_ext = args.array<double>("ebqe_bc_u_ext");
2285 xt::pyarray<int> &isDiffusiveFluxBoundary_u = args.array<int>("isDiffusiveFluxBoundary_u");
2286 xt::pyarray<int> &isAdvectiveFluxBoundary_u = args.array<int>("isAdvectiveFluxBoundary_u");
2287 xt::pyarray<double> &ebqe_bc_flux_u_ext = args.array<double>("ebqe_bc_flux_u_ext");
2288 xt::pyarray<double> &ebqe_bc_advectiveFlux_u_ext = args.array<double>("ebqe_bc_advectiveFlux_u_ext");
2289 xt::pyarray<int> &csrColumnOffsets_eb_u_u = args.array<int>("csrColumnOffsets_eb_u_u");
2290 xt::pyarray<double> &ebqe_penalty_ext = args.array<double>("ebqe_penalty_ext");
2291 const bool embeddedBoundary = args.scalar<int>("embeddedBoundary");
2292 const double embeddedBoundary_penalty = args.scalar<double>("embeddedBoundary_penalty");
2293 const double embeddedBoundary_ghost_penalty = args.scalar<double>("embeddedBoundary_ghost_penalty");
2294 xt::pyarray<double> &embeddedBoundary_sdf_nodes = args.array<double>("embeddedBoundary_sdf_nodes");
2295 xt::pyarray<double> &embeddedBoundary_sdf_q = args.array<double>("embeddedBoundary_sdf_q");
2296 xt::pyarray<double> &embeddedBoundary_normal_q = args.array<double>("embeddedBoundary_normal_q");
2297 xt::pyarray<double> &embeddedBoundary_u_q = args.array<double>("embeddedBoundary_u_q");
2298 const bool immersedBoundary = args.scalar<int>("immersedBoundary");
2299 const double immersedBoundary_penalty = args.scalar<double>("immersedBoundary_penalty");
2300 const double immersedSCIFEM_switch = args.scalar<double>("immersedSCIFEM_switch");
2301 const double immersedSCIFEM_penalty = args.scalar<double>("immersedSCIFEM_penalty");
2302 const bool PG = args.scalar<int>("PG");
2303 xt::pyarray<double> &immersedBoundary_sdf_nodes = args.array<double>("immersedBoundary_sdf_nodes");
2304 xt::pyarray<double> &immersedBoundary_sdf_q = args.array<double>("immersedBoundary_sdf_q");
2305 xt::pyarray<double> &immersedBoundary_normal_q = args.array<double>("immersedBoundary_normal_q");
2306 xt::pyarray<double> &immersedBoundary_u_q = args.array<double>("immersedBoundary_u_q");
2307 xt::pyarray<double> &immersedBoundary_fluxJump_q = args.array<double>("immersedBoundary_fluxJump_q");
2308 xt::pyarray<double> &immersedBoundary_fluxJumpVector_q = args.array<double>("immersedBoundary_fluxJumpVector_q");
2309 xt::pyarray<double> &immersedBoundary_solutionJump_nodes = args.array<double>("immersedBoundary_solutionJump_nodes");
2310 xt::pyarray<double> &isActiveDOF = args.array<double>("isActiveDOF");
2311 const double eb_adjoint_sigma = args.scalar<double>("eb_adjoint_sigma");
2312 xt::pyarray<double> &x_ref = args.array<double>("x_ref");
2313 xt::pyarray<double> &xB_ref = args.array<double>("xB_ref");
2314 xt::pyarray<int> &elementBoundariesArray = args.array<int>("elementBoundariesArray");
2315 const int nElementBoundaries_owned = args.scalar<int>("nElementBoundaries_owned");
2316 xt::pyarray<double> &elementBoundaryDiameter = args.array<double>("elementBoundaryDiameter");
2317 xt::pyarray<double> &nodeDiametersArray = args.array<double>("nodeDiametersArray");
2318 const double mua = args.scalar<double>("mua");
2319 const double mub = args.scalar<double>("mub");
2320 const double jf = args.scalar<double>("jf");
2321 const bool recomputeIFEMGeometry = args.scalar<int>("recomputeIFEMGeometry");
2322 ensureIFEMCacheSized(nElements_global); // also (re)asserts useExact = true on (re)allocation
2323 if (recomputeIFEMGeometry)
2324 ifemGeometryGeneration++;
2325 //
2326 // loop over elements to compute volume integrals and load them into the element Jacobians and global Jacobian
2327 //
2328 // std::cout << "We are computing the Jacobian...\n" << std::endl;
2329 for (int eN = 0; eN < nElements_global; eN++)
2330 {
2331 // std::cout << "########################\n element: " << eN << " \n########################" << std::endl;
2332 // double elementJacobian_u_u.data()[nDOF_test_element*nDOF_trial_element],element_u[nDOF_trial_element];
2333 auto elementJacobian_u_u = xt::pyarray<double>::from_shape({nDOF_test_element * nDOF_trial_element});
2334 auto element_u = xt::pyarray<double>::from_shape({nDOF_trial_element});
2335 for (int j = 0; j < nDOF_trial_element; j++)
2336 {
2337 int eN_j = eN * nDOF_trial_element + j;
2338 element_u.data()[j] = u_dof.data()[u_l2g.data()[eN_j]];
2339 }
2340 double element_phi_s[nDOF_trial_element];
2341 for (int j = 0; j < nDOF_trial_element; j++)
2342 {
2343 int eN_j = eN * nDOF_trial_element + j;
2344 element_phi_s[j] = embeddedBoundary_sdf_nodes.data()[u_l2g.data()[eN_j]];
2345 }
2346 double element_phi_f[nDOF_trial_element];
2347 for (int j = 0; j < nDOF_trial_element; j++)
2348 {
2349 int eN_j = eN * nDOF_trial_element + j;
2350 element_phi_f[j] = immersedBoundary_sdf_nodes.data()[u_l2g.data()[eN_j]];
2351 }
2352 double element_nodes[nDOF_trial_element * 3];
2353 for (int i = 0; i < nDOF_trial_element; i++)
2354 {
2355 int eN_i = eN * nDOF_trial_element + i;
2356 for (int I = 0; I < 3; I++)
2357 // element_nodes[i * 3 + I] = mesh_dof.data()[mesh_l2g.data()[eN_i] * 3 + I];
2358 element_nodes[i * 3 + I] = mesh_dof.data()[u_l2g.data()[eN_i] * 3 + I];
2359 // std::cout << "element_nodes[" << i << "] = (" << element_nodes[i * 3 + 0] << ", " << element_nodes[i * 3 + 1] << ", " << element_nodes[i * 3 + 2] << ")\n";
2360 } // i
2361 if (gf_s_interior_gen[eN] != ifemGeometryGeneration)
2362 {
2363 gf_s_interior_icase[eN] = gf_s_cache[eN].calculate(element_phi_s, element_nodes, x_ref.data(), false);
2364 gf_s_interior_gen[eN] = ifemGeometryGeneration;
2365 }
2366 int icase_s = gf_s_interior_icase[eN];
2367 if (gf_f_interior_gen[eN] != ifemGeometryGeneration)
2368 {
2369 gf_f_interior_icase[eN] = gf_f_cache[eN].calculate(element_phi_f, element_nodes, x_ref.data(), mua, mub, jf, false, false);
2370 gf_f_interior_gen[eN] = ifemGeometryGeneration;
2371 }
2372 int icase_f = gf_f_interior_icase[eN];
2373 calculateElementJacobian(icase_f,
2374 mesh_trial_ref,
2375 mesh_grad_trial_ref,
2376 mesh_dof,
2377 mesh_l2g,
2378 x_ref,
2379 dV_ref,
2380 u_trial_ref,
2381 u_grad_trial_ref,
2382 u_test_ref,
2383 u_grad_test_ref,
2384 elementDiameter,
2385 elementBoundaryDiameter,
2386 nodeDiametersArray,
2387 cfl,
2388 Ct_sge,
2389 sc_uref,
2390 sc_alpha,
2391 useMetrics,
2392 mesh_trial_trace_ref,
2393 mesh_grad_trial_trace_ref,
2394 dS_ref,
2395 u_trial_trace_ref,
2396 u_grad_trial_trace_ref,
2397 u_test_trace_ref,
2398 u_grad_test_trace_ref,
2399 normal_ref,
2400 boundaryJac_ref,
2401 nElements_global,
2402 nElementBoundaries_owned,
2403 u_l2g,
2404 u_dof,
2405 sd_rowptr,
2406 sd_colind,
2407 q_a,
2408 q_v,
2409 q_r,
2410 lag_shockCapturing,
2411 shockCapturingDiffusion,
2412 q_numDiff_u,
2413 q_numDiff_u_last,
2414 elementJacobian_u_u,
2415 element_u,
2416 eN,
2417 embeddedBoundary,
2418 embeddedBoundary_penalty,
2419 embeddedBoundary_normal_q,
2420 embeddedBoundary_u_q,
2421 immersedBoundary,
2422 immersedBoundary_penalty,
2423 immersedBoundary_sdf_q,
2424 immersedBoundary_normal_q,
2425 immersedBoundary_u_q,
2426 immersedBoundary_fluxJump_q,
2427 immersedBoundary_fluxJumpVector_q,
2428 immersedBoundary_solutionJump_nodes,
2429 element_phi_f,
2430 mua,
2431 mub,
2432 PG);
2433 //
2434 // load into element Jacobian into global Jacobian
2435 //
2436 for (int i = 0; i < nDOF_test_element; i++)
2437 {
2438 int eN_i = eN * nDOF_test_element + i;
2439 for (int j = 0; j < nDOF_trial_element; j++)
2440 {
2441 int eN_i_j = eN_i * nDOF_trial_element + j;
2442 globalJacobian.data()[csrRowIndeces_u_u.data()[eN_i] + csrColumnOffsets_u_u.data()[eN_i_j]] += elementJacobian_u_u.data()[i * nDOF_trial_element + j];
2443 // std::cout << "globalJacobian[" << eN_i << "," << eN * nDOF_trial_element + j << "] += " << elementJacobian_u_u.data()[i * nDOF_trial_element + j] << std::endl;
2444 } // j
2445 } // i
2446 std::cout << std::endl;
2447 } // elements
2448 for (std::set<int>::iterator it = cutfem_boundaries.begin(); it != cutfem_boundaries.end(); ++it)
2449 {
2450 std::map<int, double> Dw_Dn_jump;
2451 std::map<std::pair<int, int>, int> u_u_nz;
2452 double gamma_cutfem = embeddedBoundary_ghost_penalty, h_cutfem = elementBoundaryDiameter.data()[*it];
2453 int eN_nDOF_trial_element = elementBoundaryElementsArray.data()[(*it) * 2 + 0] * nDOF_trial_element;
2454 for (int kb = 0; kb < nQuadraturePoints_elementBoundary; kb++)
2455 {
2456 double Dp_Dn_jump = 0.0, Du_Dn_jump = 0.0, Dv_Dn_jump = 0.0, dS;
2457 for (int eN_side = 0; eN_side < 2; eN_side++)
2458 {
2459 int ebN = *it,
2460 eN = elementBoundaryElementsArray.data()[ebN * 2 + eN_side];
2461 for (int i = 0; i < nDOF_test_element; i++)
2462 Dw_Dn_jump[u_l2g.data()[eN * nDOF_test_element + i]] = 0.0;
2463 }
2464 for (int eN_side = 0; eN_side < 2; eN_side++)
2465 {
2466 int ebN = *it,
2467 eN = elementBoundaryElementsArray.data()[ebN * 2 + eN_side],
2468 ebN_local = elementBoundaryLocalElementBoundariesArray.data()[ebN * 2 + eN_side],
2469 eN_nDOF_trial_element = eN * nDOF_trial_element,
2470 ebN_local_kb = ebN_local * nQuadraturePoints_elementBoundary + kb,
2471 ebN_local_kb_nSpace = ebN_local_kb * nSpace;
2472 double
2473 u_int = 0.0,
2474 grad_u_int[nSpace] = {0., 0.},
2475 jac_int[nSpace * nSpace],
2476 jacDet_int,
2477 jacInv_int[nSpace * nSpace],
2478 boundaryJac[nSpace * (nSpace - 1)],
2479 metricTensor[(nSpace - 1) * (nSpace - 1)],
2480 metricTensorDetSqrt,
2481 u_test_dS[nDOF_test_element],
2482 u_grad_trial_trace[nDOF_trial_element * nSpace],
2483 u_grad_test_dS[nDOF_trial_element * nSpace],
2484 normal[2], x_int, y_int, z_int, xt_int, yt_int, zt_int, integralScaling,
2485 G[nSpace * nSpace], G_dd_G, tr_G, h_phi, h_penalty, penalty;
2486 // compute information about mapping from reference element to physical element
2487 ck.calculateMapping_elementBoundary(eN,
2488 ebN_local,
2489 kb,
2490 ebN_local_kb,
2491 mesh_dof.data(),
2492 mesh_l2g.data(),
2493 mesh_trial_trace_ref.data(),
2494 mesh_grad_trial_trace_ref.data(),
2495 boundaryJac_ref.data(),
2496 jac_int,
2497 jacDet_int,
2498 jacInv_int,
2499 boundaryJac,
2500 metricTensor,
2501 metricTensorDetSqrt,
2502 normal_ref.data(),
2503 normal,
2504 x_int, y_int, z_int);
2505 dS = metricTensorDetSqrt * dS_ref.data()[kb];
2506 // compute shape and solution information
2507 // shape
2508 // std::cout << "Calculating gradTrialFromRef from calculateJacobian() 1" << std::endl;
2509 ck.gradTrialFromRef(&u_grad_trial_trace_ref.data()[ebN_local_kb_nSpace * nDOF_trial_element], jacInv_int, u_grad_trial_trace);
2510 for (int i = 0; i < nDOF_test_element; i++)
2511 {
2512 int eN_i = eN * nDOF_test_element + i;
2513 for (int I = 0; I < nSpace; I++)
2514 Dw_Dn_jump[u_l2g.data()[eN_i]] += u_grad_trial_trace[i * nSpace + I] * normal[I];
2515 }
2516 } // eN_side
2517 for (int eN_side = 0; eN_side < 2; eN_side++)
2518 {
2519 int ebN = *it,
2520 eN = elementBoundaryElementsArray.data()[ebN * 2 + eN_side];
2521 for (int i = 0; i < nDOF_test_element; i++)
2522 {
2523 int eN_i = eN * nDOF_test_element + i;
2524 for (int eN_side2 = 0; eN_side2 < 2; eN_side2++)
2525 {
2526 int eN2 = elementBoundaryElementsArray.data()[ebN * 2 + eN_side2];
2527 for (int j = 0; j < nDOF_test_element; j++)
2528 {
2529 int eN_i_j = eN_i * nDOF_test_element + j;
2530 int eN2_j = eN2 * nDOF_test_element + j;
2531 int ebN_i_j = ebN * 4 * nDOF_test_X_trial_element +
2532 eN_side * 2 * nDOF_test_X_trial_element +
2533 eN_side2 * nDOF_test_X_trial_element +
2534 i * nDOF_trial_element +
2535 j;
2536 std::pair<int, int> ij = std::make_pair(u_l2g.data()[eN_i], u_l2g.data()[eN2_j]);
2537 if (u_u_nz.count(ij))
2538 {
2539 assert(u_u_nz[ij] == csrRowIndeces_u_u.data()[eN_i] + csrColumnOffsets_eb_u_u.data()[ebN_i_j]);
2540 }
2541 else
2542 u_u_nz[ij] = csrRowIndeces_u_u.data()[eN_i] + csrColumnOffsets_eb_u_u.data()[ebN_i_j];
2543 }
2544 }
2545 }
2546 }
2547 for (std::map<int, double>::iterator wi_it = Dw_Dn_jump.begin(); wi_it != Dw_Dn_jump.end(); ++wi_it)
2548 for (std::map<int, double>::iterator wj_it = Dw_Dn_jump.begin(); wj_it != Dw_Dn_jump.end(); ++wj_it)
2549 {
2550 int i_global = wi_it->first,
2551 j_global = wj_it->first;
2552 double Dw_Dn_jump_i = wi_it->second,
2553 Dw_Dn_jump_j = wj_it->second;
2554 std::pair<int, int> ij = std::make_pair(i_global, j_global);
2555 globalJacobian.data()[u_u_nz.at(ij)] += gamma_cutfem * h_cutfem * Dw_Dn_jump_j * Dw_Dn_jump_i * dS;
2556 } // i,j
2557 } // kb
2558 } // cutfem element boundaries
2559 // Per-face reference-ELEMENT quadrature points, taken from proteus' own trace tables.
2560 // Simplex::calculate(..., isBoundary=true) reads its xi_r argument as reference *element*
2561 // coordinates (it forms x = node0 + Jac_0*xi_r), but xB_ref holds reference *boundary*
2562 // points (t,0,0) -- passing those directly puts every face's points on the local
2563 // node0->node1 edge, i.e. correct only for local face 2. mesh_trial_trace_ref holds the P1
2564 // mesh shape functions at (face, quadrature point) and for a triangle those barycentric
2565 // functions ARE the reference coordinates (phi_0=1-xi-eta, phi_1=xi, phi_2=eta), so this
2566 // is exactly consistent with calculateMapping_elementBoundary's own convention.
2567 double xB_ref_faces[nDOF_mesh_trial_element * nQuadraturePoints_elementBoundary * 3];
2568 for (int f = 0; f < nDOF_mesh_trial_element; f++)
2569 for (int kb = 0; kb < nQuadraturePoints_elementBoundary; kb++)
2570 {
2571 const int fkb = f * nQuadraturePoints_elementBoundary + kb;
2572 const double *phi_m = &mesh_trial_trace_ref.data()[fkb * nDOF_mesh_trial_element];
2573 xB_ref_faces[fkb * 3 + 0] = phi_m[1];
2574 xB_ref_faces[fkb * 3 + 1] = phi_m[2];
2575 xB_ref_faces[fkb * 3 + 2] = 0.0;
2576 }
2577 // Jacobian of the SCIFEM (eq. 3.7) consistency terms assembled in calculateResidual --
2578 // same construction (orientation map + edge Heaviside blend of region-wise products),
2579 // differentiated w.r.t. each side's u_dof. Everything is linear in u_dof, so
2580 // d(avg_flux_a)/d(u_dof[s,j]) = 0.5*fta_m[s][j] and d(jump_ua)/d(u_dof[s,j]) =
2581 // (+/-)va_m[s][j] -- the same per-dof quantities used for the test index i.
2582 for (std::set<int>::iterator it = ifem_boundaries.begin(); it != ifem_boundaries.end(); ++it)
2583 {
2584 {
2585 std::map<std::pair<int, int>, int> u_u_nz;
2586 int ebN0 = *it;
2587 for (int eN_side = 0; eN_side < 2; eN_side++)
2588 {
2589 int eN = elementBoundaryElementsArray.data()[ebN0 * 2 + eN_side];
2590 for (int i = 0; i < nDOF_test_element; i++)
2591 {
2592 int eN_i = eN * nDOF_test_element + i;
2593 for (int eN_side2 = 0; eN_side2 < 2; eN_side2++)
2594 {
2595 int eN2 = elementBoundaryElementsArray.data()[ebN0 * 2 + eN_side2];
2596 for (int j = 0; j < nDOF_test_element; j++)
2597 {
2598 int eN2_j = eN2 * nDOF_test_element + j;
2599 int ebN_i_j = ebN0 * 4 * nDOF_test_X_trial_element +
2600 eN_side * 2 * nDOF_test_X_trial_element +
2601 eN_side2 * nDOF_test_X_trial_element +
2602 i * nDOF_trial_element + j;
2603 std::pair<int, int> ij = std::make_pair(u_l2g.data()[eN_i], u_l2g.data()[eN2_j]);
2604 if (!u_u_nz.count(ij))
2605 u_u_nz[ij] = csrRowIndeces_u_u.data()[eN_i] + csrColumnOffsets_eb_u_u.data()[ebN_i_j];
2606 }
2607 }
2608 }
2609 }
2610
2611 int ebN = *it;
2612 int eN_s[2], ebN_local_s[2];
2613 for (int s = 0; s < 2; s++)
2614 {
2615 eN_s[s] = elementBoundaryElementsArray.data()[ebN * 2 + s];
2616 ebN_local_s[s] = elementBoundaryLocalElementBoundariesArray.data()[ebN * 2 + s];
2617 }
2618
2619 // --- Orientation map ------------------------------------------------------
2620 // The two elements sharing an interior face do NOT necessarily parametrize it in
2621 // the same direction, so quadrature index kb can denote DIFFERENT physical points
2622 // on the two sides (measured: 4 of 14 cut edges at test=8.0/refinement=1, with
2623 // exactly the mirror offsets |1-2t_k|*L of the Gauss rule). Pairing the sides by
2624 // raw kb then differences u at two different points, so [[u]] and {{beta grad u}}
2625 // pick up u's variation ALONG the edge instead of the genuine inter-element jump.
2626 // Build an explicit map instead: kmap[s][k] is side s's index for the physical
2627 // point that side 0 calls k. Matched by position, so no assumption about the
2628 // rule being symmetric or ascending.
2629 double xq[2][nQuadraturePoints_elementBoundary][3];
2630 for (int s = 0; s < 2; s++)
2631 for (int kb = 0; kb < nQuadraturePoints_elementBoundary; kb++)
2632 {
2633 double jac_t[nSpace * nSpace], jacDet_t, jacInv_t[nSpace * nSpace],
2634 bJac_t[nSpace * (nSpace - 1)], mT_t[(nSpace - 1) * (nSpace - 1)],
2635 mTDS_t, nrm_t[2], xt_ = 0.0, yt_ = 0.0, zt_ = 0.0;
2636 ck.calculateMapping_elementBoundary(eN_s[s], ebN_local_s[s], kb,
2637 ebN_local_s[s] * nQuadraturePoints_elementBoundary + kb,
2638 mesh_dof.data(), mesh_l2g.data(),
2639 mesh_trial_trace_ref.data(), mesh_grad_trial_trace_ref.data(),
2640 boundaryJac_ref.data(), jac_t, jacDet_t, jacInv_t,
2641 bJac_t, mT_t, mTDS_t, normal_ref.data(), nrm_t,
2642 xt_, yt_, zt_);
2643 xq[s][kb][0] = xt_; xq[s][kb][1] = yt_; xq[s][kb][2] = zt_;
2644 }
2645 // edge diameter: the length scale for the pairing tolerance just below,
2646 // and for the interior-penalty scaling further down (same measure the
2647 // cutfem ghost penalty uses); note gamma/h here for a VALUE jump, vs
2648 // gamma*h there for a gradient jump.
2649 const double h_edge = elementBoundaryDiameter.data()[ebN];
2650 // Match over nSpace components only: a component the mapping does not
2651 // define contributes nothing but noise to the distance, and the paired
2652 // points agree EXACTLY in the components it does define, so the winning
2653 // distance is 0 and the runner-up is a full quadrature spacing away.
2654 // Verify that outcome rather than trust it -- a mis-pick here silently
2655 // differences u at two different points, and assert() is compiled out by
2656 // the -DNDEBUG that reaches these translation units.
2657 int kmap[2][nQuadraturePoints_elementBoundary];
2658 bool kmap_used[nQuadraturePoints_elementBoundary] = {};
2659 for (int k = 0; k < nQuadraturePoints_elementBoundary; k++)
2660 {
2661 kmap[0][k] = k;
2662 int best = -1;
2663 double bestd = 1.0e300, nextd = 1.0e300;
2664 for (int j = 0; j < nQuadraturePoints_elementBoundary; j++)
2665 {
2666 double d = 0.0;
2667 for (int I = 0; I < nSpace; I++)
2668 d += std::fabs(xq[1][j][I] - xq[0][k][I]);
2669 if (d < bestd) { nextd = bestd; bestd = d; best = j; }
2670 else if (d < nextd) { nextd = d; }
2671 }
2672 if (best < 0 || kmap_used[best] || bestd >= 1.0e-8 * h_edge || nextd <= 2.0 * bestd)
2673 {
2674 std::cerr << "ADR immersedSCIFEM: quadrature points on face " << ebN
2675 << " between elements " << eN_s[0] << " and " << eN_s[1]
2676 << " do not pair up (k=" << k << ", best=" << best
2677 << ", d=" << bestd << ", runner-up=" << nextd
2678 << ", h=" << h_edge << "); the mesh is non-conforming or the"
2679 << " boundary mapping is inconsistent." << std::endl;
2680 throw std::runtime_error("ADR immersedSCIFEM: face quadrature pairing failed");
2681 }
2682 kmap_used[best] = true;
2683 kmap[1][k] = best;
2684 }
2685
2686 // --- Edge Heaviside moment fit -------------------------------------------
2687 // The facet integrand is DISCONTINUOUS at the interface crossing theta, so a fixed
2688 // whole-edge Gauss rule cannot integrate it directly. Weight the two one-sided
2689 // integrands by an edge-restricted moment-fit Heaviside instead: exact for any
2690 // polynomial integrand of degree <= nP_edge, hence exact for P1 (degree <= 1) and
2691 // P2 (degree <= 3) alike -- no linearity assumption, unlike a split-quadrature
2692 // scheme. ONE fit per edge, in side 0's parametrization: the two sides' fits are
2693 // mirror images, moment-equivalent but not pointwise equal, so they must not be
2694 // mixed. t below is therefore always measured along side 0's edge direction.
2695 double edge_nodes[2][3], phi_edge[2];
2696 {
2697 const int n0l = (ebN_local_s[0] + 1) % 3, n1l = (ebN_local_s[0] + 2) % 3;
2698 const int ln[2] = {n0l, n1l};
2699 for (int e = 0; e < 2; e++)
2700 {
2701 int eN_i = eN_s[0] * nDOF_trial_element + ln[e];
2702 phi_edge[e] = immersedBoundary_sdf_nodes.data()[u_l2g.data()[eN_i]];
2703 for (int I = 0; I < 3; I++)
2704 edge_nodes[e][I] = mesh_dof.data()[u_l2g.data()[eN_i] * 3 + I];
2705 }
2706 }
2707 double C_H_edge[nP_edge + 1];
2708 equivalent_polynomials::calculate_edge_H<nP_edge>(phi_edge[0], phi_edge[1], C_H_edge);
2709 double edge_vec[3];
2710 double edge_len2 = 0.0;
2711 for (int I = 0; I < 3; I++)
2712 {
2713 edge_vec[I] = edge_nodes[1][I] - edge_nodes[0][I];
2714 edge_len2 += edge_vec[I] * edge_vec[I];
2715 }
2716
2717 for (int k = 0; k < nQuadraturePoints_elementBoundary; k++)
2718 {
2719 double dS = 0.0;
2720 std::map<int, double> va_m[2], vb_m[2], fta_m[2], ftb_m[2];
2721 for (int eN_side = 0; eN_side < 2; eN_side++)
2722 {
2723 int eN = eN_s[eN_side],
2724 ebN_local = ebN_local_s[eN_side],
2725 kb = kmap[eN_side][k],
2726 eN_nDOF_trial_element = eN * nDOF_trial_element,
2727 ebN_local_kb = ebN_local * nQuadraturePoints_elementBoundary + kb;
2728 double jac_int[nSpace * nSpace], jacDet_int, jacInv_int[nSpace * nSpace],
2729 boundaryJac[nSpace * (nSpace - 1)],
2730 metricTensor[(nSpace - 1) * (nSpace - 1)], metricTensorDetSqrt,
2731 normal[2], x_int, y_int, z_int;
2732 ck.calculateMapping_elementBoundary(eN, ebN_local, kb, ebN_local_kb,
2733 mesh_dof.data(), mesh_l2g.data(),
2734 mesh_trial_trace_ref.data(), mesh_grad_trial_trace_ref.data(),
2735 boundaryJac_ref.data(), jac_int, jacDet_int, jacInv_int,
2736 boundaryJac, metricTensor, metricTensorDetSqrt,
2737 normal_ref.data(), normal, x_int, y_int, z_int);
2738 if (eN_side == 0)
2739 dS = metricTensorDetSqrt * dS_ref.data()[kb];
2740 double sign_ne = (eN_side == 0) ? 1.0 : -1.0;
2741
2742 double element_phi_f[nDOF_trial_element], element_nodes[nDOF_trial_element * 3];
2743 for (int i = 0; i < nDOF_trial_element; i++)
2744 {
2745 int eN_i = eN * nDOF_trial_element + i;
2746 element_phi_f[i] = immersedBoundary_sdf_nodes.data()[u_l2g.data()[eN_i]];
2747 for (int I = 0; I < 3; I++)
2748 element_nodes[i * 3 + I] = mesh_dof.data()[u_l2g.data()[eN_i] * 3 + I];
2749 }
2750 if (gf_f_boundary_gen[eN] != ifemGeometryGeneration)
2751 {
2752 gf_f_boundary_icase[eN] = gf_f_cache[eN].calculate(element_phi_f, element_nodes, xB_ref_faces, mua, mub, jf, true, false);
2753 gf_f_boundary_gen[eN] = ifemGeometryGeneration;
2754 }
2755 GfType &gf_f = gf_f_cache[eN];
2756 gf_f.set_boundary_quad(ebN_local_kb);
2757
2758 for (int i = 0; i < nDOF_test_element; i++)
2759 {
2760 int dof_i = u_l2g.data()[eN_nDOF_trial_element + i];
2761 va_m[eN_side][dof_i] = gf_f.VA(i);
2762 vb_m[eN_side][dof_i] = gf_f.VB(i);
2763 double gax = gf_f.VA_x(i), gay = gf_f.VA_y(i),
2764 gbx = gf_f.VB_x(i), gby = gf_f.VB_y(i);
2765 fta_m[eN_side][dof_i] = sign_ne * mua * (gax * normal[0] + gay * normal[1]);
2766 ftb_m[eN_side][dof_i] = sign_ne * mub * (gbx * normal[0] + gby * normal[1]);
2767 }
2768 } // eN_side
2769
2770 // t of this physical point along side 0's edge direction, then H_hat(t)
2771 double proj = 0.0;
2772 for (int I = 0; I < 3; I++)
2773 proj += (xq[0][k][I] - edge_nodes[0][I]) * edge_vec[I];
2774 double t_edge = (edge_len2 > 1.0e-24) ? (proj / edge_len2) : 0.0;
2775 double H_e = equivalent_polynomials::evaluate_edge_poly<nP_edge>(C_H_edge, t_edge);
2776 double ImH_e = 1.0 - H_e;
2777
2778 for (int side_i = 0; side_i < 2; side_i++)
2779 {
2780 double vi_sign = (side_i == 0) ? 1.0 : -1.0;
2781 for (std::map<int, double>::iterator wi = va_m[side_i].begin(); wi != va_m[side_i].end(); ++wi)
2782 {
2783 int dof_i = wi->first;
2784 double jva_i = vi_sign * va_m[side_i][dof_i], jvb_i = vi_sign * vb_m[side_i][dof_i];
2785 double fta_i = fta_m[side_i][dof_i], ftb_i = ftb_m[side_i][dof_i];
2786 for (int side_j = 0; side_j < 2; side_j++)
2787 {
2788 double vj_sign = (side_j == 0) ? 1.0 : -1.0;
2789 for (std::map<int, double>::iterator wj = va_m[side_j].begin(); wj != va_m[side_j].end(); ++wj)
2790 {
2791 int dof_j = wj->first;
2792 double jva_j = vj_sign * va_m[side_j][dof_j], jvb_j = vj_sign * vb_m[side_j][dof_j];
2793 double fta_j = fta_m[side_j][dof_j], ftb_j = ftb_m[side_j][dof_j];
2794 double dPa = 0.5 * fta_j * jva_i + 0.5 * fta_i * jva_j;
2795 double dPb = 0.5 * ftb_j * jvb_i + 0.5 * ftb_i * jvb_j;
2796 std::pair<int, int> ij = std::make_pair(dof_i, dof_j);
2797 globalJacobian.data()[u_u_nz.at(ij)] -=
2798 immersedSCIFEM_switch * (ImH_e * dPa + H_e * dPb) * dS;
2799 // derivative of the interior-penalty term: d[[u]]/du_j = [[v_j]],
2800 // so the contribution is the (symmetric, PSD) gram term below.
2801 globalJacobian.data()[u_u_nz.at(ij)] +=
2802 (immersedSCIFEM_penalty / h_edge) *
2803 (ImH_e * jva_j * jva_i + H_e * jvb_j * jvb_i) * dS;
2804 }
2805 }
2806 }
2807 }
2808 } // k
2809 }
2810 } // ifem element boundaries
2811 //
2812 // loop over exterior element boundaries to compute the surface integrals and load them into the global Jacobian
2813 //
2814 for (int ebNE = 0; ebNE < nExteriorElementBoundaries_global; ebNE++)
2815 {
2816 int ebN = exteriorElementBoundariesArray.data()[ebNE];
2817 int eN = elementBoundaryElementsArray.data()[ebN * 2 + 0],
2818 ebN_local = elementBoundaryLocalElementBoundariesArray.data()[ebN * 2 + 0],
2819 eN_nDOF_trial_element = eN * nDOF_trial_element;
2820 for (int kb = 0; kb < nQuadraturePoints_elementBoundary; kb++)
2821 {
2822 int ebNE_kb = ebNE * nQuadraturePoints_elementBoundary + kb,
2823 ebNE_kb_nSpace = ebNE_kb * nSpace,
2824 ebN_local_kb = ebN_local * nQuadraturePoints_elementBoundary + kb,
2825 ebN_local_kb_nSpace = ebN_local_kb * nSpace;
2826
2827 double u_ext = 0.0,
2828 grad_u_ext[nSpace],
2829 m_ext = 0.0,
2830 dm_ext = 0.0,
2831 *a_ext,
2832 f_ext[nSpace],
2833 df_ext[nSpace],
2834 r_ext = 0.0,
2835 dflux_u_u_ext = 0.0,
2836 bc_u_ext = 0.0,
2837 // bc_grad_u_ext[nSpace],
2838 bc_m_ext = 0.0,
2839 bc_dm_ext = 0.0,
2840 bc_f_ext[nSpace],
2841 bc_df_ext[nSpace],
2842 fluxJacobian_u_u[nDOF_trial_element],
2843 jac_ext[nSpace * nSpace],
2844 jacDet_ext,
2845 jacInv_ext[nSpace * nSpace],
2846 boundaryJac[nSpace * (nSpace - 1)],
2847 metricTensor[(nSpace - 1) * (nSpace - 1)],
2848 metricTensorDetSqrt,
2849 dS,
2850 u_test_dS[nDOF_test_element],
2851 u_grad_trial_trace[nDOF_trial_element * nSpace],
2852 u_grad_test_dS[nDOF_trial_element * nSpace],
2853 normal[nSpace], x_ext, y_ext, z_ext, xt_ext, yt_ext, zt_ext, integralScaling,
2854 //
2855 G[nSpace * nSpace], G_dd_G, tr_G;
2856 ck.calculateMapping_elementBoundary(eN,
2857 ebN_local,
2858 kb,
2859 ebN_local_kb,
2860 mesh_dof.data(),
2861 mesh_l2g.data(),
2862 mesh_trial_trace_ref.data(),
2863 mesh_grad_trial_trace_ref.data(),
2864 boundaryJac_ref.data(),
2865 jac_ext,
2866 jacDet_ext,
2867 jacInv_ext,
2868 boundaryJac,
2869 metricTensor,
2870 metricTensorDetSqrt,
2871 normal_ref.data(),
2872 normal,
2873 x_ext, y_ext, z_ext);
2874 dS = metricTensorDetSqrt * dS_ref.data()[kb];
2875 ck.calculateG(jacInv_ext, G, G_dd_G, tr_G);
2876 // compute shape and solution information
2877 // shape
2878 // std::cout << "Calculating gradTrialFromRef from calculateJacobian() 3" << std::endl;
2879 ck.gradTrialFromRef(&u_grad_trial_trace_ref.data()[ebN_local_kb_nSpace * nDOF_trial_element], jacInv_ext, u_grad_trial_trace);
2880 // solution and gradients
2881 ck.valFromDOF(u_dof.data(), &u_l2g.data()[eN_nDOF_trial_element], &u_trial_trace_ref.data()[ebN_local_kb * nDOF_test_element], u_ext);
2882 ck.gradFromDOF(u_dof.data(), &u_l2g.data()[eN_nDOF_trial_element], u_grad_trial_trace, grad_u_ext);
2883 // precalculate test function products with integration weights
2884 for (int j = 0; j < nDOF_trial_element; j++)
2885 {
2886 u_test_dS[j] = u_test_trace_ref.data()[ebN_local_kb * nDOF_test_element + j] * dS;
2887 for (int I = 0; I < nSpace; I++)
2888 u_grad_test_dS[j * nSpace + I] = u_grad_trial_trace[j * nSpace + I] * dS; // cek hack, using trial
2889 }
2890 //
2891 // load the boundary values
2892 //
2893 bc_u_ext = isDOFBoundary_u.data()[ebNE_kb] * ebqe_bc_u_ext.data()[ebNE_kb] + (1 - isDOFBoundary_u.data()[ebNE_kb]) * u_ext;
2894 a_ext = &ebqe_a.data()[ebNE_kb * sd_rowptr.data()[nSpace]];
2895 for (int I = 0; I < nSpace; I++)
2896 {
2897 df_ext[I] = ebqe_v.data()[ebNE_kb * nSpace + I];
2898 bc_df_ext[I] = ebqe_v.data()[ebNE_kb * nSpace + I];
2899 }
2900 //
2901 // calculate the numerical fluxes
2902 //
2903 exteriorNumericalAdvectiveFluxDerivative(isDOFBoundary_u.data()[ebNE_kb],
2904 isAdvectiveFluxBoundary_u.data()[ebNE_kb],
2905 normal,
2906 df_ext,
2907 dflux_u_u_ext);
2908 //
2909 // calculate the flux jacobian
2910 //
2911 for (int j = 0; j < nDOF_trial_element; j++)
2912 {
2913 // int ebNE_kb_j = ebNE_kb*nDOF_trial_element+j;
2914 int j_nSpace = j * nSpace, ebN_local_kb_j = ebN_local_kb * nDOF_trial_element + j;
2915 fluxJacobian_u_u[j] = ExteriorNumericalDiffusiveFluxJacobian(sd_rowptr.data(),
2916 sd_colind.data(),
2917 isDOFBoundary_u.data()[ebNE_kb],
2918 isDiffusiveFluxBoundary_u.data()[ebNE_kb],
2919 normal,
2920 a_ext,
2921 u_trial_trace_ref.data()[ebN_local_kb_j],
2922 &u_grad_trial_trace[j_nSpace],
2923 ebqe_penalty_ext.data()[ebNE_kb]) +
2924 ck.ExteriorNumericalAdvectiveFluxJacobian(dflux_u_u_ext, u_trial_trace_ref.data()[ebN_local_kb_j]);
2925 } // j
2926 //
2927 // update the global Jacobian from the flux Jacobian
2928 //
2929 for (int i = 0; i < nDOF_test_element; i++)
2930 {
2931 int eN_i = eN * nDOF_test_element + i;
2932 int i_nSpace = i * nSpace;
2933 for (int j = 0; j < nDOF_trial_element; j++)
2934 {
2935 int ebN_i_j = ebN * 4 * nDOF_test_X_trial_element + i * nDOF_trial_element + j;
2936 int ebN_local_kb_j = ebN_local_kb * nDOF_trial_element + j;
2937
2938 globalJacobian.data()[csrRowIndeces_u_u.data()[eN_i] + csrColumnOffsets_eb_u_u.data()[ebN_i_j]] += fluxJacobian_u_u[j] * u_test_dS[i] +
2939 ck.ExteriorElementBoundaryDiffusionAdjointJacobian(isDOFBoundary_u.data()[ebNE_kb],
2940 isDiffusiveFluxBoundary_u.data()[ebNE_kb],
2941 eb_adjoint_sigma,
2942 u_trial_trace_ref.data()[ebN_local_kb_j],
2943 normal,
2944 sd_rowptr.data(),
2945 sd_colind.data(),
2946 a_ext,
2947 &u_grad_test_dS[i_nSpace]);
2948 } // j
2949 } // i
2950 } // kb
2951 } // ebNE
2952 } // computeJacobian
2953 }; // cADR
2954
2955 inline cADR_base *newADR(int nSpaceIn,
2956 int nQuadraturePoints_elementIn,
2957 int nDOF_mesh_trial_elementIn,
2958 int nDOF_trial_elementIn,
2959 int nDOF_test_elementIn,
2960 int nQuadraturePoints_elementBoundaryIn,
2961 int CompKernelFlag)
2962 {
2963 if (nSpaceIn == 2)
2964 return proteus::chooseAndAllocateDiscretization2D<cADR_base, cADR, CompKernel>(nSpaceIn,
2965 nQuadraturePoints_elementIn,
2966 nDOF_mesh_trial_elementIn,
2967 nDOF_trial_elementIn,
2968 nDOF_test_elementIn,
2969 nQuadraturePoints_elementBoundaryIn,
2970 CompKernelFlag);
2971 else
2972 return proteus::chooseAndAllocateDiscretization<cADR_base, cADR, CompKernel>(nSpaceIn,
2973 nQuadraturePoints_elementIn,
2974 nDOF_mesh_trial_elementIn,
2975 nDOF_trial_elementIn,
2976 nDOF_test_elementIn,
2977 nQuadraturePoints_elementBoundaryIn,
2978 CompKernelFlag);
2979 }
2980} // proteus
2981#endif
n
Int n
Definition Headers.h:28
sol
Double ** sol
Definition Headers.h:40
r
Double r
Definition Headers.h:83
s
Double s
Definition Headers.h:84
u
Double u
Definition Headers.h:89
num
Int num
Definition Headers.h:32
z
Double * z
Definition Headers.h:49
v
Double v
Definition Headers.h:95
equivalent_polynomials::GeneralizedFunctions_mix
Definition equivalent_polynomials.h:1285
equivalent_polynomials::GeneralizedFunctions_mix::VA_y
double VA_y(int i)
Definition equivalent_polynomials.h:1374
equivalent_polynomials::GeneralizedFunctions_mix::VA_x
double VA_x(int i)
Definition equivalent_polynomials.h:1367
equivalent_polynomials::GeneralizedFunctions_mix::set_quad
void set_quad(unsigned int q)
Definition equivalent_polynomials.h:1325
equivalent_polynomials::GeneralizedFunctions_mix::VA
double VA(int i)
Definition equivalent_polynomials.h:1360
equivalent_polynomials::GeneralizedFunctions_mix::H
double H(double eps, double phi)
Definition equivalent_polynomials.h:1337
equivalent_polynomials::GeneralizedFunctions_mix::exact
Simplex< nSpace, nP_ifem, nP, nQ, nEBQ, useIfemBasis > exact
Definition equivalent_polynomials.h:1288
equivalent_polynomials::GeneralizedFunctions_mix::ImH
double ImH(double eps, double phi)
Definition equivalent_polynomials.h:1345
equivalent_polynomials::GeneralizedFunctions_mix::D
double D(double eps, double phi)
Definition equivalent_polynomials.h:1353
equivalent_polynomials::GeneralizedFunctions_mix::VB_y
double VB_y(int i)
Definition equivalent_polynomials.h:1402
equivalent_polynomials::GeneralizedFunctions_mix::VB
double VB(int i)
Definition equivalent_polynomials.h:1388
equivalent_polynomials::GeneralizedFunctions_mix::VB_x
double VB_x(int i)
Definition equivalent_polynomials.h:1395
equivalent_polynomials::GeneralizedFunctions_mix::set_boundary_quad
void set_boundary_quad(unsigned int ebq)
Definition equivalent_polynomials.h:1331
equivalent_polynomials::GeneralizedFunctions_mix::get_normal
double * get_normal()
Definition equivalent_polynomials.h:1317
proteus::cADR_base
Definition ADR.h:24
proteus::cADR_base::calculateResidual
virtual void calculateResidual(arguments_dict &args)=0
proteus::cADR_base::calculateJacobian
virtual void calculateJacobian(arguments_dict &args)=0
proteus::cADR_base::~cADR_base
virtual ~cADR_base()
Definition ADR.h:26
proteus::cADR::updateImmersedBoundaryTerms
void updateImmersedBoundaryTerms(const double immersedBoundary_penalty, const double dV, double *immersedBoundary_normal, double x, double y, double z, const double u_s, const double u, const double grad_u[nSpace], const double a, double &r, double &dr, double &ham, double *dham, double *f, double *df, const double fluxJump, const double *fluxJumpVector, const double D_f)
Definition ADR.h:335
proteus::cADR::exteriorNumericalDiffusiveFlux
void exteriorNumericalDiffusiveFlux(int *rowptr, int *colind, const int &isDOFBoundary, const int &isDiffusiveFluxBoundary, const double n[nSpace], double *bc_a, const double &bc_u, const double &bc_flux, double *a, const double grad_potential[nSpace], const double &u, const double &penalty, double &flux)
Definition ADR.h:91
proteus::cADR::gf_s_interior_gen
std::vector< int > gf_s_interior_gen
Definition ADR.h:67
proteus::cADR::calculateSubgridError_tau
void calculateSubgridError_tau(const double &Ct_sge, const double G[nSpace *nSpace], const double &A0, const double Ai[nSpace], double &tau_v, double &cfl)
Definition ADR.h:185
proteus::cADR::gf_s_interior_icase
std::vector< int > gf_s_interior_icase
Definition ADR.h:68
proteus::cADR::nP_edge
static const int nP_edge
Definition ADR.h:54
proteus::cADR::gf_s_cache
std::vector< GfType > gf_s_cache
Definition ADR.h:66
proteus::cADR::ifem_boundaries
std::set< int > ifem_boundaries
Definition ADR.h:41
proteus::cADR::ck
CompKernelType ck
Definition ADR.h:45
proteus::cADR::calculateSubgridError_tau
void calculateSubgridError_tau(const double &elementDiameter, const double &dmt, const double dH[nSpace], double &cfl, double &tau)
Definition ADR.h:166
proteus::cADR::nDOF_test_X_trial_element
const int nDOF_test_X_trial_element
Definition ADR.h:44
proteus::cADR::ifem_boundary_elements
std::set< int > ifem_boundary_elements
Definition ADR.h:41
proteus::cADR::ensureIFEMCacheSized
void ensureIFEMCacheSized(int nElements_global)
Definition ADR.h:73
proteus::cADR::exteriorNumericalAdvectiveFluxDerivative
void exteriorNumericalAdvectiveFluxDerivative(const int &isDOFBoundary_u, const int &isFluxBoundary_u, const double n[nSpace], const double velocity[nSpace], double &dflux)
Definition ADR.h:262
proteus::cADR::cutfem_boundary_elements
std::set< int > cutfem_boundary_elements
Definition ADR.h:42
proteus::cADR::ifemGeometryGeneration
int ifemGeometryGeneration
Definition ADR.h:69
proteus::cADR::elementIsActive
std::valarray< bool > elementIsActive
Definition ADR.h:43
proteus::cADR::calculateElementResidual
void calculateElementResidual(int icase_f, xt::pyarray< double > &mesh_trial_ref, xt::pyarray< double > &mesh_grad_trial_ref, xt::pyarray< double > &mesh_dof, xt::pyarray< int > &mesh_l2g, xt::pyarray< double > &x_ref, xt::pyarray< double > &dV_ref, xt::pyarray< double > &u_trial_ref, xt::pyarray< double > &u_grad_trial_ref, xt::pyarray< double > &u_test_ref, xt::pyarray< double > &u_grad_test_ref, xt::pyarray< double > &elementDiameter, xt::pyarray< double > &elementBoundaryDiameter, xt::pyarray< double > &nodeDiametersArray, xt::pyarray< double > &cfl, double Ct_sge, double sc_uref, double sc_alpha, double useMetrics, xt::pyarray< double > &mesh_trial_trace_ref, xt::pyarray< double > &mesh_grad_trial_trace_ref, xt::pyarray< double > &dS_ref, xt::pyarray< double > &u_trial_trace_ref, xt::pyarray< double > &u_grad_trial_trace_ref, xt::pyarray< double > &u_test_trace_ref, xt::pyarray< double > &u_grad_test_trace_ref, xt::pyarray< double > &normal_ref, xt::pyarray< double > &boundaryJac_ref, int nElements_global, int nElementBoundaries_owned, xt::pyarray< int > &u_l2g, xt::pyarray< double > &u_dof, xt::pyarray< int > &sd_rowptr, xt::pyarray< int > &sd_colind, xt::pyarray< double > &q_a, xt::pyarray< double > &q_v, xt::pyarray< double > &q_r, int lag_shockCapturingDiffusion, double shockCapturingDiffusion, xt::pyarray< double > &q_numDiff_u, xt::pyarray< double > &q_numDiff_u_last, int offset_u, int stride_u, xt::pyarray< double > &elementResidual_u, int nExteriorElementBoundaries_global, xt::pyarray< int > &exteriorElementBoundariesArray, xt::pyarray< int > &elementBoundariesArray, xt::pyarray< int > &elementBoundaryElementsArray, xt::pyarray< int > &elementBoundaryLocalElementBoundariesArray, xt::pyarray< double > &element_u, int eN, const bool embeddedBoundary, const double embeddedBoundary_penalty, xt::pyarray< double > &embeddedBoundary_normal_q, xt::pyarray< double > &embeddedBoundary_u_q, const bool immersedBoundary, const double immersedBoundary_penalty, xt::pyarray< double > &immersedBoundary_sdf_q, xt::pyarray< double > &immersedBoundary_normal_q, xt::pyarray< double > &immersedBoundary_u_q, xt::pyarray< double > &immersedBoundary_fluxJump_q, xt::pyarray< double > &immersedBoundary_fluxJumpVector_q, xt::pyarray< double > &immersedBoundary_solutionJump_nodes, double *element_phi_f, bool &element_active, std::valarray< bool > &elementIsActive, double *JA, double *JB, double &L2_error, double &Linfty_error, double mua, double mub, xt::pyarray< double > &q_u_exact_inner, xt::pyarray< double > &q_u_exact_outer, const bool PG)
Definition ADR.h:394
proteus::cADR::calculateElementJacobian
void calculateElementJacobian(int icase_f, xt::pyarray< double > &mesh_trial_ref, xt::pyarray< double > &mesh_grad_trial_ref, xt::pyarray< double > &mesh_dof, xt::pyarray< int > &mesh_l2g, xt::pyarray< double > &x_ref, xt::pyarray< double > &dV_ref, xt::pyarray< double > &u_trial_ref, xt::pyarray< double > &u_grad_trial_ref, xt::pyarray< double > &u_test_ref, xt::pyarray< double > &u_grad_test_ref, xt::pyarray< double > &elementDiameter, xt::pyarray< double > &elementBoundaryDiameter, xt::pyarray< double > &nodeDiametersArray, xt::pyarray< double > &cfl, double Ct_sge, double sc_uref, double sc_alpha, double useMetrics, xt::pyarray< double > &mesh_trial_trace_ref, xt::pyarray< double > &mesh_grad_trial_trace_ref, xt::pyarray< double > &dS_ref, xt::pyarray< double > &u_trial_trace_ref, xt::pyarray< double > &u_grad_trial_trace_ref, xt::pyarray< double > &u_test_trace_ref, xt::pyarray< double > &u_grad_test_trace_ref, xt::pyarray< double > &normal_ref, xt::pyarray< double > &boundaryJac_ref, int nElements_global, int nElementBoundaries_owned, xt::pyarray< int > &u_l2g, xt::pyarray< double > &u_dof, xt::pyarray< int > &sd_rowptr, xt::pyarray< int > &sd_colind, xt::pyarray< double > &q_a, xt::pyarray< double > &q_v, xt::pyarray< double > &q_r, int lag_shockCapturing, double shockCapturingDiffusion, xt::pyarray< double > &q_numDiff_u, xt::pyarray< double > &q_numDiff_u_last, xt::pyarray< double > &elementJacobian_u_u, xt::pyarray< double > &element_u, int eN, const bool embeddedBoundary, const double embeddedBoundary_penalty, xt::pyarray< double > &embeddedBoundary_normal_q, xt::pyarray< double > &embeddedBoundary_u_q, const bool immersedBoundary, const double immersedBoundary_penalty, xt::pyarray< double > &immersedBoundary_sdf_q, xt::pyarray< double > &immersedBoundary_normal_q, xt::pyarray< double > &immersedBoundary_u_q, xt::pyarray< double > &immersedBoundary_fluxJump_q, xt::pyarray< double > &immersedBoundary_fluxJumpVector_q, xt::pyarray< double > &immersedBoundary_solutionJump_nodes, double *element_phi_f, double mua, double mub, const bool PG)
Definition ADR.h:1817
proteus::cADR::gf_f_boundary_gen
std::vector< int > gf_f_boundary_gen
Definition ADR.h:67
proteus::cADR::gf_f_interior_icase
std::vector< int > gf_f_interior_icase
Definition ADR.h:68
proteus::cADR::gf_f_boundary_icase
std::vector< int > gf_f_boundary_icase
Definition ADR.h:68
proteus::cADR::calculateNumericalDiffusion
void calculateNumericalDiffusion(const double &shockCapturingDiffusion, const double &elementDiameter, const double &strong_residual, const double grad_u[nSpace], double &numDiff)
Definition ADR.h:199
proteus::cADR::updateEmbeddedBoundaryTerms
void updateEmbeddedBoundaryTerms(const double embeddedBoundary_penalty, const double dV, double *embeddedBoundary_normal, const double u_s, const double u, const double grad_u[nSpace], const double a, double &r, double &dr, double &ham, double *dham, double *f, double *df, const double D_s)
Definition ADR.h:299
proteus::cADR::cutfem_boundaries
std::set< int > cutfem_boundaries
Definition ADR.h:42
proteus::cADR::ExteriorNumericalDiffusiveFluxJacobian
double ExteriorNumericalDiffusiveFluxJacobian(int *rowptr, int *colind, const int &isDOFBoundary, const int &isDiffusiveFluxBoundary, const double n[nSpace], double *a, const double &v, const double grad_v[nSpace], const double &penalty)
Definition ADR.h:136
proteus::cADR::gf_f_interior_gen
std::vector< int > gf_f_interior_gen
Definition ADR.h:67
proteus::cADR::calculateResidual
void calculateResidual(arguments_dict &args)
Definition ADR.h:924
proteus::cADR::GfType
GeneralizedFunctions< nSpace, nDOF_trial_element, 4, nQuadraturePoints_element, nDOF_mesh_trial_element *nQuadraturePoints_elementBoundary > GfType
Definition ADR.h:50
proteus::cADR::calculateJacobian
void calculateJacobian(arguments_dict &args)
Definition ADR.h:2236
proteus::cADR::exteriorNumericalAdvectiveFlux
void exteriorNumericalAdvectiveFlux(const int &isDOFBoundary_u, const int &isFluxBoundary_u, const double n[nSpace], const double &bc_u, const double &bc_flux_u, const double &u, const double velocity[nSpace], double &flux)
Definition ADR.h:218
proteus::cADR::gf_f_cache
std::vector< GfType > gf_f_cache
Definition ADR.h:66
equivalent_polynomials.h
df
double df(double C, double b, double a, int q, int r)
Definition analyticalSolutions.c:2209
sign
#define sign(x, y)
Definition jf.h:44
w
#define w(x)
Definition jf.h:22
equivalent_polynomials::evaluate_edge_poly
double evaluate_edge_poly(const double C[nP+1], double t)
Definition equivalent_polynomials_utils.h:151
equivalent_polynomials::calculate_edge_H
void calculate_edge_H(double phi0, double phi1, double C_H[nP+1])
Definition equivalent_polynomials_utils.h:132
proteus
Definition ADR.h:19
proteus::GeneralizedFunctions
equivalent_polynomials::GeneralizedFunctions_mix< nSpace, nP_ifem, nP, nQ, nEBQ, true > GeneralizedFunctions
Definition ADR.h:21
proteus::newADR
cADR_base * newADR(int nSpaceIn, int nQuadraturePoints_elementIn, int nDOF_mesh_trial_elementIn, int nDOF_trial_elementIn, int nDOF_test_elementIn, int nQuadraturePoints_elementBoundaryIn, int CompKernelFlag)
Definition ADR.h:2955
proteus::chooseAndAllocateDiscretization
Model_Base * chooseAndAllocateDiscretization(int nSpaceIn, int nQuadraturePoints_elementIn, int nDOF_mesh_trial_elementIn, int nDOF_trial_elementIn, int nDOF_test_elementIn, int nDOF_v_trial_elementIn, int nDOF_v_test_elementIn, int nQuadraturePoints_elementBoundaryIn, int CompKernelFlag)
Definition MixedModelFactory.h:48
proteus::chooseAndAllocateDiscretization2D
Model_Base * chooseAndAllocateDiscretization2D(int nSpaceIn, int nQuadraturePoints_elementIn, int nDOF_mesh_trial_elementIn, int nDOF_trial_elementIn, int nDOF_test_elementIn, int nDOF_v_trial_elementIn, int nDOF_v_test_elementIn, int nQuadraturePoints_elementBoundaryIn, int CompKernelFlag)
Definition MixedModelFactory.h:309
proteus::f
double f(const double &g, const double &h, const double &hZ)
Definition SW2DCV.h:58
proteus::arguments_dict
Definition ArgumentsDict.h:70
proteus::arguments_dict::scalar
T & scalar(const std::string &key)
proteus::arguments_dict::array
xt::pyarray< T > & array(const std::string &key)
Generated on for proteus by doxygen 1.18.0