proteus  1.8.1
C/C++/Fortran libraries
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 <map>
7 #include <valarray>
8 #include "CompKernel.h"
9 #include "ModelFactory.h"
10 #include "equivalent_polynomials.h"
11 #include "mprans/ArgumentsDict.h"
12 #include "xtensor-python/pyarray.hpp"
13 
14 namespace py = pybind11;
15 
16 namespace proteus
17 {
18  template<int nSpace, int nP, int nQ, int nEBQ>
19  using GeneralizedFunctions = equivalent_polynomials::GeneralizedFunctions_mix<nSpace, nP, nQ, nEBQ>;
20 
21  class cADR_base
22  {
23  public:
24  virtual ~cADR_base(){}
25  virtual void calculateResidual(arguments_dict& args) = 0;
26  virtual void calculateJacobian(arguments_dict& args) = 0;
27  };
28 
29  template<class CompKernelType,
30  int nSpace,
31  int nQuadraturePoints_element,
32  int nDOF_mesh_trial_element,
33  int nDOF_trial_element,
34  int nDOF_test_element,
35  int nQuadraturePoints_elementBoundary
36  >
37  class cADR: public cADR_base
38  {
39  public:
40  std::set<int> ifem_boundaries, ifem_boundary_elements,
41  cutfem_boundaries, cutfem_boundary_elements;
42  std::valarray<bool> elementIsActive;
43  const int nDOF_test_X_trial_element;
44  CompKernelType ck;
45  GeneralizedFunctions<nSpace,3,nQuadraturePoints_element,nQuadraturePoints_elementBoundary> gf;
46  GeneralizedFunctions<nSpace,3,nQuadraturePoints_element,nQuadraturePoints_elementBoundary> gf_s;
47  cADR(): nDOF_test_X_trial_element(nDOF_test_element*nDOF_trial_element), ck()
48  {}
49  /* inline
50  void evaluateCoefficients()
51  {
52  }
53  */
54 
55  inline
56  void exteriorNumericalDiffusiveFlux(int* rowptr,
57  int* colind,
58  const int& isDOFBoundary,
59  const int& isDiffusiveFluxBoundary,
60  const double n[nSpace],
61  double* bc_a,
62  const double& bc_u,
63  const double& bc_flux,
64  double* a,
65  const double grad_potential[nSpace],
66  const double& u,
67  const double& penalty,
68  double& flux)
69  {
70  double diffusiveVelocityComponent_I;
71  double penaltyFlux;
72  double max_a;
73  if (isDiffusiveFluxBoundary == 1)
74  {
75  flux = bc_flux;
76  }
77  else if (isDOFBoundary == 1)
78  {
79  flux = 0.0;
80  max_a = 0.0;
81  for (int I = 0; I < nSpace; I++)
82  {
83  diffusiveVelocityComponent_I = 0.0;
84  for (int m = rowptr[I]; m < rowptr[I+1]; m++)
85  {
86  diffusiveVelocityComponent_I -= a[m] * grad_potential[colind[m]];
87  max_a = fmax(max_a, a[m]);
88  }
89  flux += diffusiveVelocityComponent_I * n[I];
90  }
91  penaltyFlux = max_a * penalty * (u - bc_u);
92  flux += penaltyFlux;
93  }
94  else
95  {
96  std::cerr << "warning, diffusion term with no boundary condition set, setting diffusive flux to 0.0" << std::endl;
97  flux = 0.0;
98  }
99  }
100 
101  inline
102  double ExteriorNumericalDiffusiveFluxJacobian(int* rowptr,
103  int* colind,
104  const int& isDOFBoundary,
105  const int& isDiffusiveFluxBoundary,
106  const double n[nSpace],
107  double* a,
108  const double& v,
109  const double grad_v[nSpace],
110  const double& penalty)
111  {
112  double dvel_I;
113  double tmp = 0.0;
114  double max_a = 0.0;
115  if ((isDiffusiveFluxBoundary == 0) && (isDOFBoundary == 1))
116  {
117  for (int I = 0; I < nSpace; I++)
118  {
119  dvel_I = 0.0;
120  for (int m = rowptr[I]; m < rowptr[I + 1]; m++)
121  {
122  dvel_I -= a[m] * grad_v[colind[m]];
123  max_a = fmax(max_a, a[m]);
124  }
125  tmp += dvel_I * n[I];
126  }
127  tmp += max_a * penalty * v;
128  }
129  return tmp;
130  }
131 
132  inline
133  void calculateSubgridError_tau(const double& elementDiameter,
134  const double& dmt,
135  const double dH[nSpace],
136  double& cfl,
137  double& tau)
138  {
139  double h;
140  double nrm_v;
141  double oneByAbsdt;
142  h = elementDiameter;
143  nrm_v = 0.0;
144  for(int I = 0; I < nSpace; I++)
145  nrm_v += dH[I] * dH[I];
146  nrm_v = sqrt(nrm_v);
147  cfl = nrm_v / h;
148  oneByAbsdt = fabs(dmt);
149  tau = 1.0 / (2.0 * nrm_v / h + oneByAbsdt + 1.0e-8);
150  }
151 
152  inline
153  void calculateSubgridError_tau(const double& Ct_sge,
154  const double G[nSpace*nSpace],
155  const double& A0,
156  const double Ai[nSpace],
157  double& tau_v,
158  double& cfl)
159  {
160  double v_d_Gv = 0.0;
161  for (int I = 0; I < nSpace; I++)
162  for (int J = 0; J < nSpace; J++)
163  v_d_Gv += Ai[I] * G[I * nSpace + J] * Ai[J];
164  tau_v = 1.0 / sqrt(Ct_sge * A0 * A0 + v_d_Gv + 1.0e-8);
165  }
166 
167  inline
168  void calculateNumericalDiffusion(const double& shockCapturingDiffusion,
169  const double& elementDiameter,
170  const double& strong_residual,
171  const double grad_u[nSpace],
172  double& numDiff)
173  {
174  double h;
175  double num;
176  double den;
177  double n_grad_u;
178  h = elementDiameter;
179  n_grad_u = 0.0;
180  for (int I = 0; I < nSpace; I++)
181  n_grad_u += grad_u[I] * grad_u[I];
182  num = shockCapturingDiffusion * 0.5 * h * fabs(strong_residual);
183  den = sqrt(n_grad_u) + 1.0e-8;
184  numDiff = num / den;
185  }
186 
187  inline
188  void exteriorNumericalAdvectiveFlux(const int& isDOFBoundary_u,
189  const int& isFluxBoundary_u,
190  const double n[nSpace],
191  const double& bc_u,
192  const double& bc_flux_u,
193  const double& u,
194  const double velocity[nSpace],
195  double& flux)
196  {
197 
198  double flow = 0.0;
199  for (int I = 0; I < nSpace; I++)
200  flow += n[I] * velocity[I];
201  if (isDOFBoundary_u == 1)
202  {
203  if (flow >= 0.0)
204  {
205  flux = u * flow;
206  //flux = flow;
207  }
208  else
209  {
210  flux = bc_u * flow;
211  //flux = flow;
212  }
213  }
214  else if (isFluxBoundary_u == 1)
215  {
216  flux = bc_flux_u;
217  }
218  else
219  {
220  if (flow >= 0.0)
221  {
222  flux = u * flow;
223  }
224  else
225  {
226  flux = 0.0;
227  }
228 
229  }
230  //flux = flow;
231  }
232 
233  inline
234  void exteriorNumericalAdvectiveFluxDerivative(const int& isDOFBoundary_u,
235  const int& isFluxBoundary_u,
236  const double n[nSpace],
237  const double velocity[nSpace],
238  double& dflux)
239  {
240  double flow = 0.0;
241  for (int I = 0; I < nSpace; I++)
242  {
243  flow += n[I] * velocity[I];
244  }
245  //double flow=n[0]*velocity[0]+n[1]*velocity[1]+n[2]*velocity[2];
246  dflux = 0.0;//default to no flux
247  if (isDOFBoundary_u == 1)
248  {
249  if (flow >= 0.0)
250  {
251  dflux = flow;
252  }
253  else
254  {
255  dflux = 0.0;
256  }
257  }
258  else if (isFluxBoundary_u == 1)
259  {
260  dflux = 0.0;
261  }
262  else
263  {
264  if (flow >= 0.0)
265  {
266  dflux = flow;
267  }
268  }
269  }
270 
271  inline void updateEmbeddedBoundaryTerms(const double embeddedBoundary_penalty,
272  const double dV,
273  double* embeddedBoundary_normal,
274  const double u_s,
275  const double u,
276  const double grad_u[nSpace],
277  const double a,
278  double& r,
279  double& dr,
280  double& ham,
281  double* dham,
282  double* f,
283  double* df)
284  {
285  double outward_normal[2]={0.0,0.0};
286  for (int I=0;I<nSpace;I++)
287  outward_normal[I] = -embeddedBoundary_normal[I];
288 
289  double D_s = gf_s.D(0.,0.);
290 
291  //diffusive flux
292  ham -= D_s * a * (outward_normal[0] * grad_u[0] + outward_normal[1] * grad_u[1]);
293  dham[0] -= D_s * a * outward_normal[0];
294  dham[1] -= D_s * a * outward_normal[1];
295 
296  //Nitsche Dirichlet penalty
297  r += D_s*a*embeddedBoundary_penalty * (u - u_s);
298  dr += D_s*a*embeddedBoundary_penalty;
299 
300  //Nitsche adjoint consistency
301  f[0] += D_s * a * outward_normal[0] * (u - u_s);
302  f[1] += D_s * a * outward_normal[1] * (u - u_s);
303  df[0] += D_s * a * outward_normal[0];
304  df[1] += D_s * a * outward_normal[1];
305  }
306 
307  inline void calculateElementResidual(//element
308  xt::pyarray<double>& mesh_trial_ref,
309  xt::pyarray<double>& mesh_grad_trial_ref,
310  xt::pyarray<double>& mesh_dof,
311  xt::pyarray<int>& mesh_l2g,
312  xt::pyarray<double>& x_ref,
313  xt::pyarray<double>& dV_ref,
314  xt::pyarray<double>& u_trial_ref,
315  xt::pyarray<double>& u_grad_trial_ref,
316  xt::pyarray<double>& u_test_ref,
317  xt::pyarray<double>& u_grad_test_ref,
318  xt::pyarray<double>& elementDiameter,
319  xt::pyarray<double>& elementBoundaryDiameter,
320  xt::pyarray<double>& nodeDiametersArray,
321  xt::pyarray<double>& cfl,
322  double Ct_sge,
323  double sc_uref,
324  double sc_alpha,
325  double useMetrics,
326  //element boundary
327  xt::pyarray<double>& mesh_trial_trace_ref,
328  xt::pyarray<double>& mesh_grad_trial_trace_ref,
329  xt::pyarray<double>& dS_ref,
330  xt::pyarray<double>& u_trial_trace_ref,
331  xt::pyarray<double>& u_grad_trial_trace_ref,
332  xt::pyarray<double>& u_test_trace_ref,
333  xt::pyarray<double>& u_grad_test_trace_ref,
334  xt::pyarray<double>& normal_ref,
335  xt::pyarray<double>& boundaryJac_ref,
336  //physics
337  int nElements_global,
338  int nElementBoundaries_owned,
339  xt::pyarray<int>& u_l2g,
340  xt::pyarray<double>& u_dof,
341  xt::pyarray<int>& sd_rowptr,
342  xt::pyarray<int>& sd_colind,
343  xt::pyarray<double>& q_a,
344  xt::pyarray<double>& q_v,
345  xt::pyarray<double>& q_r,
346  int lag_shockCapturingDiffusion,
347  double shockCapturingDiffusion,
348  xt::pyarray<double>& q_numDiff_u,
349  xt::pyarray<double>& q_numDiff_u_last,
350  int offset_u,
351  int stride_u,
352  xt::pyarray<double>& elementResidual_u,
353  int nExteriorElementBoundaries_global,
354  xt::pyarray<int>& exteriorElementBoundariesArray,
355  xt::pyarray<int>& elementBoundariesArray,
356  xt::pyarray<int>& elementBoundaryElementsArray,
357  xt::pyarray<int>& elementBoundaryLocalElementBoundariesArray,
358  xt::pyarray<double>& element_u,
359  int eN,
360  const bool embeddedBoundary,
361  const double embeddedBoundary_penalty,
362  xt::pyarray<double>& embeddedBoundary_normal_q,
363  xt::pyarray<double>& embeddedBoundary_u_q,
364  bool& element_active,
365  std::valarray<bool>& elementIsActive)
366  {
367  for (int i = 0; i < nDOF_test_element; i++)
368  {
369  elementResidual_u.data()[i] = 0.0;
370  }
371  //loop over quadrature points and compute integrands
372  for (int k = 0; k < nQuadraturePoints_element; k++)
373  {
374  gf_s.set_quad(k);
375  //compute indeces and declare local storage
376  int eN_k = eN * nQuadraturePoints_element + k;
377  int eN_k_3d = eN_k*3;
378  double h_phi;
379  double u = 0.0;
380  double grad_u[nSpace];
381  double m = 0.0;
382  double dm=0.0;
383  double f[nSpace];
384  double df[nSpace];
385  double f_s[nSpace]={0.,0.};
386  double df_s[nSpace]={0.,0.};
387  double ham_s=0.0;
388  double dham_s[nSpace]={0.,0.};
389  double m_t = 0.0;
390  double dm_t = 0.0;
391  double pdeResidual_u = 0.0;
392  double Lstar_u[nDOF_test_element];
393  double subgridError_u = 0.0;
394  double tau = 0.0;
395  double tau0 = 0.0;
396  double tau1 = 0.0;
397  double numDiff0 = 0.0;
398  double numDiff1 = 0.0;
399  double *a = NULL;
400  double r = 0.0;
401  double r_s = 0.0;
402  double dr_s = 0.0;
403  double jac[nSpace*nSpace];
404  double jacDet;
405  double jacInv[nSpace*nSpace];
406  double u_grad_trial[nDOF_trial_element * nSpace];
407  double u_test_dV[nDOF_trial_element];
408  double u_grad_test_dV[nDOF_test_element * nSpace];
409  double dV;
410  double x;
411  double y;
412  double z;
413  double G[nSpace * nSpace];
414  double G_dd_G;
415  double tr_G;
416  //
417  //compute solution and gradients at quadrature points
418  //
419  ck.calculateMapping_element(eN,
420  k,
421  mesh_dof.data(),
422  mesh_l2g.data(),
423  mesh_trial_ref.data(),
424  mesh_grad_trial_ref.data(),
425  jac,
426  jacDet,
427  jacInv,
428  x,
429  y,
430  z);
431  ck.calculateH_element(eN,
432  k,
433  nodeDiametersArray.data(),
434  mesh_l2g.data(),
435  mesh_trial_ref.data(),
436  h_phi);
437  //get the physical integration weight
438  dV = fabs(jacDet) * dV_ref.data()[k];
439  //get the metric tensor and friends
440  ck.calculateG(jacInv, G, G_dd_G, tr_G);
441  //get the trial function gradients
442  ck.gradTrialFromRef(&u_grad_trial_ref.data()[k * nDOF_trial_element * nSpace], jacInv, u_grad_trial);
443  //get the solution
444  ck.valFromElementDOF(element_u.data(), &u_trial_ref.data()[k * nDOF_trial_element], u);
445  //get the solution gradients
446  ck.gradFromElementDOF(element_u.data(), u_grad_trial, grad_u);
447  //precalculate test function products with integration weights
448  for (int j = 0; j < nDOF_trial_element; j++)
449  {
450  u_test_dV[j] = u_test_ref.data()[k * nDOF_trial_element + j] * dV;
451  for (int I = 0; I < nSpace; I++)
452  {
453  u_grad_test_dV[j * nSpace + I] = u_grad_trial[j * nSpace + I] *dV; //cek warning won't work for Petrov-Galerkin
454  }
455  }
456  //
457  //calculate pde coefficients at quadrature points
458  //
459  //evaluateCoefficients();
460  //just set from pre-evaluated quadrature point values for now
461  a = &q_a.data()[eN_k * sd_rowptr.data()[nSpace]];
462  r = q_r.data()[eN_k];
463  for (int I = 0; I < nSpace; I++)
464  {
465  f[I] = q_v.data()[eN_k * nSpace + I] * u;
466  df[I] = q_v.data()[eN_k * nSpace + I];
467  }
468  const double H_s = gf_s.H(0.,0.);
469  const double D_s = gf_s.D(0.,0.);
470  if ( H_s != 0.0 || D_s != 0.0)
471  {
472  element_active=true;
473  elementIsActive[eN]=true;
474  }
475  if(embeddedBoundary)
476  {
477  double level_set_normal[nSpace];
478  double sign=0.0;
479  double norm_exact=0.0,norm_cut=0.0;
480  for (int I=0;I<nSpace;I++)
481  {
482  sign += embeddedBoundary_normal_q.data()[eN_k_3d+I]*gf_s.get_normal()[I];
483  level_set_normal[I] = gf_s.get_normal()[I];
484  norm_cut += level_set_normal[I]*level_set_normal[I];
485  norm_exact += embeddedBoundary_normal_q.data()[eN_k_3d+I]*embeddedBoundary_normal_q.data()[eN_k_3d+I];
486  }
487  assert(std::fabs(1.0-norm_cut) < 1.0e-8);
488  assert(std::fabs(1.0-norm_exact) < 1.0e-8);
489  if (sign < 0.0)
490  for (int I=0;I<nSpace;I++)
491  level_set_normal[I]*=-1.0;
492  updateEmbeddedBoundaryTerms(embeddedBoundary_penalty/h_phi,//penalty,
493  dV,
494  level_set_normal,
495  embeddedBoundary_u_q.data()[eN_k],
496  u,
497  grad_u,
498  a[0],//assume scalar diffusion for now
499  r_s,
500  dr_s,
501  ham_s,
502  dham_s,
503  f_s,
504  df_s);
505  }
506  //
507  //moving mesh
508  //
509  /* double mesh_velocity[3]; */
510  /* mesh_velocity[0] = xt; */
511  /* mesh_velocity[1] = yt; */
512  /* mesh_velocity[2] = zt; */
513  /* for (int I=0;I<nSpace;I++) */
514  /* { */
515  /* f[I] -= MOVING_DOMAIN*m*mesh_velocity[I]; */
516  /* df[I] -= MOVING_DOMAIN*dm*mesh_velocity[I]; */
517  /* } */
518  //
519  //calculate time derivative at quadrature points
520  //
521  /* ck.bdf(alphaBDF, */
522  /* q_m_betaBDF.data()[eN_k], */
523  /* m, */
524  /* dm, */
525  /* m_t, */
526  /* dm_t); */
527  //
528  //calculate subgrid error (strong residual and adjoint)
529  //
530  //calculate strong residual
531  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);
532  //calculate adjoint
533  for (int i = 0; i < nDOF_test_element; i++)
534  {
535  // register int eN_k_i_nSpace = (eN_k*nDOF_trial_element+i)*nSpace;
536  // Lstar_u[i] = ck.Advection_adjoint(df,&u_grad_test_dV.data()[eN_k_i_nSpace]);
537  int i_nSpace = i * nSpace;
538  Lstar_u[i] = ck.Advection_adjoint(df, &u_grad_test_dV[i_nSpace]);
539  }
540  //calculate tau and tau*Res
541  calculateSubgridError_tau(elementDiameter.data()[eN], dm_t, df, cfl.data()[eN_k], tau0);
542  calculateSubgridError_tau(Ct_sge,
543  G,
544  dm_t,
545  df,
546  tau1,
547  cfl.data()[eN_k]);
548 
549  tau = useMetrics * tau1 + (1.0 - useMetrics) * tau0;
550 
551  subgridError_u = -tau * pdeResidual_u;
552  //
553  //calculate shock capturing diffusion
554  //
555  ck.calculateNumericalDiffusion(shockCapturingDiffusion, elementDiameter.data()[eN], pdeResidual_u, grad_u, numDiff0);
556  ck.calculateNumericalDiffusion(shockCapturingDiffusion, sc_uref, sc_alpha, G, G_dd_G, pdeResidual_u, grad_u, numDiff1);
557  q_numDiff_u.data()[eN_k] = useMetrics * numDiff1 + (1.0 - useMetrics) * numDiff0;
558  //
559  //update element residual
560  //
561  for (int i = 0; i < nDOF_test_element; i++)
562  {
563  int i_nSpace=i*nSpace;
564  elementResidual_u.data()[i] += H_s*(ck.Advection_weak(f, &u_grad_test_dV[i_nSpace]) +
565  ck.Diffusion_weak(sd_rowptr.data(), sd_colind.data(), a, grad_u, &u_grad_test_dV[i_nSpace]) +
566  ck.Reaction_weak(r, u_test_dV[i]) +
567  ck.SubgridError(subgridError_u, Lstar_u[i]) +
568  ck.NumericalDiffusion(q_numDiff_u_last.data()[eN_k], grad_u, &u_grad_test_dV[i_nSpace]));
569  if (embeddedBoundary)
570  {
571  elementResidual_u.data()[i] += (ck.Advection_weak(f_s,&u_grad_test_dV[i_nSpace]) +
572  ck.Reaction_weak(r_s,u_test_dV[i]) +
573  ck.Hamiltonian_weak(ham_s,u_test_dV[i]));
574  }
575  }
576  }
577  }
578 
579  void calculateResidual(arguments_dict& args)
580  {
581  xt::pyarray<double>& mesh_trial_ref = args.array<double>("mesh_trial_ref");
582  xt::pyarray<double>& mesh_grad_trial_ref = args.array<double>("mesh_grad_trial_ref");
583  xt::pyarray<double>& mesh_dof = args.array<double>("mesh_dof");
584  xt::pyarray<int>& mesh_l2g = args.array<int>("mesh_l2g");
585  xt::pyarray<double>& dV_ref = args.array<double>("dV_ref");
586  xt::pyarray<double>& u_trial_ref = args.array<double>("u_trial_ref");
587  xt::pyarray<double>& u_grad_trial_ref = args.array<double>("u_grad_trial_ref");
588  xt::pyarray<double>& u_test_ref = args.array<double>("u_test_ref");
589  xt::pyarray<double>& u_grad_test_ref = args.array<double>("u_grad_test_ref");
590  xt::pyarray<double>& elementDiameter = args.array<double>("elementDiameter");
591  xt::pyarray<double>& cfl = args.array<double>("cfl");
592  double Ct_sge = args.scalar<double>("Ct_sge");
593  double sc_uref = args.scalar<double>("sc_uref");
594  double sc_alpha = args.scalar<double>("sc_alpha");
595  double useMetrics = args.scalar<double>("useMetrics");
596  xt::pyarray<double>& mesh_trial_trace_ref = args.array<double>("mesh_trial_trace_ref");
597  xt::pyarray<double>& mesh_grad_trial_trace_ref = args.array<double>("mesh_grad_trial_trace_ref");
598  xt::pyarray<double>& dS_ref = args.array<double>("dS_ref");
599  xt::pyarray<double>& u_trial_trace_ref = args.array<double>("u_trial_trace_ref");
600  xt::pyarray<double>& u_grad_trial_trace_ref = args.array<double>("u_grad_trial_trace_ref");
601  xt::pyarray<double>& u_test_trace_ref = args.array<double>("u_test_trace_ref");
602  xt::pyarray<double>& u_grad_test_trace_ref = args.array<double>("u_grad_test_trace_ref");
603  xt::pyarray<double>& normal_ref = args.array<double>("normal_ref");
604  xt::pyarray<double>& boundaryJac_ref = args.array<double>("boundaryJac_ref");
605  int nElements_global = args.scalar<int>("nElements_global");
606  xt::pyarray<int>& u_l2g = args.array<int>("u_l2g");
607  xt::pyarray<double>& u_dof = args.array<double>("u_dof");
608  xt::pyarray<int>& sd_rowptr = args.array<int>("sd_rowptr");
609  xt::pyarray<int>& sd_colind = args.array<int>("sd_colind");
610  xt::pyarray<double>& q_a = args.array<double>("q_a");
611  xt::pyarray<double>& q_v = args.array<double>("q_v");
612  xt::pyarray<double>& q_r = args.array<double>("q_r");
613  int lag_shockCapturing = args.scalar<int>("lag_shockCapturing");
614  double shockCapturingDiffusion = args.scalar<double>("shockCapturingDiffusion");
615  xt::pyarray<double>& q_numDiff_u = args.array<double>("q_numDiff_u");
616  xt::pyarray<double>& q_numDiff_u_last = args.array<double>("q_numDiff_u_last");
617  int offset_u = args.scalar<int>("offset_u");
618  int stride_u = args.scalar<int>("stride_u");
619  xt::pyarray<double>& globalResidual = args.array<double>("globalResidual");
620  int nExteriorElementBoundaries_global = args.scalar<int>("nExteriorElementBoundaries_global");
621  xt::pyarray<int>& exteriorElementBoundariesArray = args.array<int>("exteriorElementBoundariesArray");
622  xt::pyarray<int>& elementBoundaryElementsArray = args.array<int>("elementBoundaryElementsArray");
623  xt::pyarray<int>& elementBoundaryLocalElementBoundariesArray = args.array<int>("elementBoundaryLocalElementBoundariesArray");
624  xt::pyarray<double>& ebqe_a = args.array<double>("ebqe_a");
625  xt::pyarray<double>& ebqe_v = args.array<double>("ebqe_v");
626  xt::pyarray<int>& isDOFBoundary_u = args.array<int>("isDOFBoundary_u");
627  xt::pyarray<double>& ebqe_bc_u_ext = args.array<double>("ebqe_bc_u_ext");
628  xt::pyarray<int>& isDiffusiveFluxBoundary_u = args.array<int>("isDiffusiveFluxBoundary_u");
629  xt::pyarray<int>& isAdvectiveFluxBoundary_u = args.array<int>("isAdvectiveFluxBoundary_u");
630  xt::pyarray<double>& ebqe_bc_flux_u_ext = args.array<double>("ebqe_bc_flux_u_ext");
631  xt::pyarray<double>& ebqe_bc_advectiveFlux_u_ext = args.array<double>("ebqe_bc_advectiveFlux_u_ext");
632  xt::pyarray<double>& ebqe_penalty_ext = args.array<double>("ebqe_penalty_ext");
633  const bool embeddedBoundary = args.scalar<int>("embeddedBoundary");
634  const double embeddedBoundary_penalty = args.scalar<double>("embeddedBoundary_penalty");
635  const double embeddedBoundary_ghost_penalty = args.scalar<double>("embedded_ghost_penalty");
636  xt::pyarray<double>& embeddedBoundary_sdf_nodes = args.array<double>("embeddedBoundary_sdf_nodes");
637  xt::pyarray<double>& embeddedBoundary_sdf_q = args.array<double>("embeddedBoundary_sdf_q");
638  xt::pyarray<double>& embeddedBoundary_normal_q = args.array<double>("embeddedBoundary_normal_q");
639  xt::pyarray<double>& embeddedBoundary_u_q = args.array<double>("embeddedBoundary_u_q");
640  xt::pyarray<double>& isActiveDOF = args.array<double>("isActiveDOF");
641  const double eb_adjoint_sigma = args.scalar<double>("eb_adjoint_sigma");
642  xt::pyarray<double>& x_ref = args.array<double>("x_ref");
643  xt::pyarray<int>& elementBoundariesArray = args.array<int>("elementBoundariesArray");
644  const int nElementBoundaries_owned = args.scalar<int>("nElementBoundaries_owned");
645  xt::pyarray<double>& elementBoundaryDiameter = args.array<double>("elementBoundaryDiameter");
646  xt::pyarray<double>& nodeDiametersArray = args.array<double>("nodeDiametersArray");
647 
648  gf.useExact = true;
649  gf_s.useExact = true;
650  ifem_boundaries.clear();
651  ifem_boundary_elements.clear();
652  cutfem_boundaries.clear();
653  cutfem_boundary_elements.clear();
654  elementIsActive.resize(nElements_global);
655 
656  //
657  //loop over elements to compute volume integrals and load them into element and global residual
658  //
659  //eN is the element index
660  //eN_k is the quadrature point index for a scalar
661  //eN_k_nSpace is the quadrature point index for a vector
662  //eN_i is the element test function index
663  //eN_j is the element trial function index
664  //eN_k_j is the quadrature point index for a trial function
665  //eN_k_i is the quadrature point index for a trial function
666  for(int eN=0;eN<nElements_global;eN++)
667  {
668  //declare local storage for element residual and initialize
669  //register double elementResidual_u[nDOF_test_element],element_u[nDOF_trial_element];
670  auto elementResidual_u = xt::pyarray<double>::from_shape({nDOF_test_element});
671  auto element_u = xt::pyarray<double>::from_shape({nDOF_trial_element});
672  bool element_active=false;
673  elementIsActive[eN]=false;
674  for (int i=0;i<nDOF_test_element;i++)
675  {
676  register int eN_i=eN*nDOF_test_element+i;
677  element_u.data()[i] = u_dof.data()[u_l2g.data()[eN_i]];
678  }//i
679  double element_phi_s[nDOF_mesh_trial_element];
680  for (int j=0;j<nDOF_mesh_trial_element;j++)
681  {
682  register int eN_j = eN*nDOF_mesh_trial_element+j;
683  element_phi_s[j] = embeddedBoundary_sdf_nodes.data()[u_l2g.data()[eN_j]];
684  }
685  double element_nodes[nDOF_mesh_trial_element*3];
686  for (int i=0;i<nDOF_mesh_trial_element;i++)
687  {
688  register int eN_i=eN*nDOF_mesh_trial_element+i;
689  for(int I=0;I<3;I++)
690  element_nodes[i*3 + I] = mesh_dof.data()[mesh_l2g.data()[eN_i]*3 + I];
691  }//i
692  int icase_s = gf_s.calculate(element_phi_s, element_nodes, x_ref.data(),false);
693  if (icase_s == 0)
694  {
695  //only works for simplices
696  for (int ebN_element=0;ebN_element < nDOF_mesh_trial_element; ebN_element++)
697  {
698  const int ebN = elementBoundariesArray.data()[eN*nDOF_mesh_trial_element+ebN_element];
699  //internal and actually a cut edge
700  //if (elementBoundaryElementsArray[ebN*2+1] != -1 && element_phi_s[(ebN_element+1)%nDOF_mesh_trial_element]*element_phi_s[(ebN_element+2)%nDOF_mesh_trial_element] <= 0.0)
701  // cutfem_boundaries.insert(ebN);
702  if (elementBoundaryElementsArray.data()[ebN*2+1] != -1 && (ebN < nElementBoundaries_owned))
703  cutfem_boundaries.insert(ebN);
704  }
705  }
706  // int icase = gf.calculate(element_phi, element_nodes, x_ref.data(), rho_1*nu_1, rho_0*nu_0,false,false);
707  // if (icase == 0)
708  // {
709  // //only works for simplices
710  // for (int ebN_element=0;ebN_element < nDOF_mesh_trial_element; ebN_element++)
711  // {
712  // const int ebN = elementBoundariesArray.data()[eN*nDOF_mesh_trial_element+ebN_element];
713  // ifem_boundaries.insert(ebN);
714  // }
715  // }
716  calculateElementResidual(mesh_trial_ref,
717  mesh_grad_trial_ref,
718  mesh_dof,
719  mesh_l2g,
720  x_ref,
721  dV_ref,
722  u_trial_ref,
723  u_grad_trial_ref,
724  u_test_ref,
725  u_grad_test_ref,
726  elementDiameter,
727  elementBoundaryDiameter,
728  nodeDiametersArray,
729  cfl,
730  Ct_sge,
731  sc_uref,
732  sc_alpha,
733  useMetrics,
734  mesh_trial_trace_ref,
735  mesh_grad_trial_trace_ref,
736  dS_ref,
737  u_trial_trace_ref,
738  u_grad_trial_trace_ref,
739  u_test_trace_ref,
740  u_grad_test_trace_ref,
741  normal_ref,
742  boundaryJac_ref,
743  nElements_global,
744  nElementBoundaries_owned,
745  u_l2g,
746  u_dof,
747  sd_rowptr,
748  sd_colind,
749  q_a,
750  q_v,
751  q_r,
752  lag_shockCapturing,
753  shockCapturingDiffusion,
754  q_numDiff_u,
755  q_numDiff_u_last,
756  offset_u,stride_u,
757  elementResidual_u,
758  nExteriorElementBoundaries_global,
759  exteriorElementBoundariesArray,
760  elementBoundariesArray,
761  elementBoundaryElementsArray,
762  elementBoundaryLocalElementBoundariesArray,
763  element_u,
764  eN,
765  embeddedBoundary,
766  embeddedBoundary_penalty,
767  embeddedBoundary_normal_q,
768  embeddedBoundary_u_q,
769  element_active,
770  elementIsActive);
771  //
772  //load element into global residual and save element residual
773  //
774  for(int i=0;i<nDOF_test_element;i++)
775  {
776  register int eN_i=eN*nDOF_test_element+i;
777 
778  globalResidual.data()[offset_u+stride_u*u_l2g.data()[eN_i]]+=elementResidual_u.data()[i];
779  if (element_active)
780  isActiveDOF.data()[offset_u+stride_u*u_l2g.data()[eN_i]] = 1.0;
781  }//i
782  }//elements
783  for (std::set<int>::iterator it=cutfem_boundaries.begin(); it!=cutfem_boundaries.end(); )
784  {
785  if(elementIsActive[elementBoundaryElementsArray[(*it)*2+0]] && elementIsActive[elementBoundaryElementsArray[(*it)*2+1]])
786  {
787  std::map<int,double> Dwp_Dn_jump, Dw_Dn_jump;
788  register double gamma_cutfem=embeddedBoundary_ghost_penalty,h_cutfem=elementBoundaryDiameter.data()[*it];
789  for (int kb=0;kb<nQuadraturePoints_elementBoundary;kb++)
790  {
791  register double Du_Dn_jump=0.0, dS;
792  for (int eN_side=0;eN_side < 2; eN_side++)
793  {
794  register int ebN = *it,
795  eN = elementBoundaryElementsArray.data()[ebN*2+eN_side];
796  for (int i=0;i<nDOF_test_element;i++)
797  {
798  Dw_Dn_jump[u_l2g.data()[eN*nDOF_test_element+i]] = 0.0;
799  }
800  }
801  for (int eN_side=0;eN_side < 2; eN_side++)
802  {
803  register int ebN = *it,
804  eN = elementBoundaryElementsArray.data()[ebN*2+eN_side],
805  ebN_local = elementBoundaryLocalElementBoundariesArray.data()[ebN*2+eN_side],
806  eN_nDOF_trial_element = eN*nDOF_trial_element,
807  ebN_local_kb = ebN_local*nQuadraturePoints_elementBoundary+kb,
808  ebN_local_kb_nSpace = ebN_local_kb*nSpace;
809  register double u_int=0.0,
810  grad_u_int[nSpace]={0.,0.},
811  jac_int[nSpace*nSpace],
812  jacDet_int,
813  jacInv_int[nSpace*nSpace],
814  boundaryJac[nSpace*(nSpace-1)],
815  metricTensor[(nSpace-1)*(nSpace-1)],
816  metricTensorDetSqrt,
817  u_test_dS[nDOF_test_element],
818  u_grad_trial_trace[nDOF_trial_element*nSpace],
819  u_grad_test_dS[nDOF_trial_element*nSpace],
820  normal[2],x_int,y_int,z_int,xt_int,yt_int,zt_int,integralScaling,
821  G[nSpace*nSpace],G_dd_G,tr_G,h_phi,h_penalty,penalty,
822  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;
823  //compute information about mapping from reference element to physical element
824  ck.calculateMapping_elementBoundary(eN,
825  ebN_local,
826  kb,
827  ebN_local_kb,
828  mesh_dof.data(),
829  mesh_l2g.data(),
830  mesh_trial_trace_ref.data(),
831  mesh_grad_trial_trace_ref.data(),
832  boundaryJac_ref.data(),
833  jac_int,
834  jacDet_int,
835  jacInv_int,
836  boundaryJac,
837  metricTensor,
838  metricTensorDetSqrt,
839  normal_ref.data(),
840  normal,
841  x_int,y_int,z_int);
842  dS = metricTensorDetSqrt*dS_ref.data()[kb];
843  //compute shape and solution information
844  //shape
845  ck.gradTrialFromRef(&u_grad_trial_trace_ref.data()[ebN_local_kb_nSpace*nDOF_trial_element],jacInv_int,u_grad_trial_trace);
846  //solution and gradients
847  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);
848  ck.gradFromDOF(u_dof.data(),&u_l2g.data()[eN_nDOF_trial_element],u_grad_trial_trace,grad_u_int);
849  for (int I=0;I<nSpace;I++)
850  {
851  Du_Dn_jump += grad_u_int[I]*normal[I];
852  }
853  for (int i=0;i<nDOF_test_element;i++)
854  {
855  for (int I=0;I<nSpace;I++)
856  Dw_Dn_jump[u_l2g.data()[eN_nDOF_trial_element+i]] += u_grad_trial_trace[i*nSpace+I]*normal[I];
857  }
858  }//eN_side
859  for (std::map<int,double>::iterator w_it=Dw_Dn_jump.begin(); w_it!=Dw_Dn_jump.end(); ++w_it)
860  {
861  int i_global = w_it->first;
862  double Dw_Dn_jump_i = w_it->second;
863  globalResidual.data()[offset_u+stride_u*i_global]+=gamma_cutfem*h_cutfem*Du_Dn_jump*Dw_Dn_jump_i*dS;
864  }//i
865  }//kb
866  ++it;
867  }
868  else
869  {
870  it = cutfem_boundaries.erase(it);
871  }
872  }//cutfem element boundaries
873  //
874  //loop over exterior element boundaries to calculate surface integrals and load into element and global residuals
875  //
876  //ebNE is the Exterior element boundary INdex
877  //ebN is the element boundary INdex
878  //eN is the element index
879  for (int ebNE = 0; ebNE < nExteriorElementBoundaries_global; ebNE++)
880  {
881  register int ebN = exteriorElementBoundariesArray.data()[ebNE],
882  eN = elementBoundaryElementsArray.data()[ebN*2+0],
883  ebN_local = elementBoundaryLocalElementBoundariesArray.data()[ebN*2+0],
884  eN_nDOF_trial_element = eN*nDOF_trial_element;
885  register double elementResidual_u[nDOF_test_element];
886  for (int i=0;i<nDOF_test_element;i++)
887  {
888  elementResidual_u[i]=0.0;
889  }
890  for (int kb=0;kb<nQuadraturePoints_elementBoundary;kb++)
891  {
892  register int ebNE_kb = ebNE*nQuadraturePoints_elementBoundary+kb,
893  ebNE_kb_nSpace = ebNE_kb*nSpace,
894  ebN_local_kb = ebN_local*nQuadraturePoints_elementBoundary+kb,
895  ebN_local_kb_nSpace = ebN_local_kb*nSpace;
896  register double u_ext=0.0,
897  grad_u_ext[nSpace],
898  m_ext=0.0,
899  dm_ext=0.0,
900  *a_ext,
901  /* *da_exxt, */
902  f_ext[nSpace],
903  df_ext[nSpace],
904  r_ext=0.0,
905  /* dr_ext=0.0, */
906  flux_diff_ext=0.0,
907  flux_advect_ext=0.0,
908  bc_u_ext=0.0,
909  //bc_grad_u_ext[nSpace],
910  bc_m_ext=0.0,
911  bc_dm_ext=0.0,
912  bc_f_ext[nSpace],
913  bc_df_ext[nSpace],
914  jac_ext[nSpace*nSpace],
915  jacDet_ext,
916  jacInv_ext[nSpace*nSpace],
917  boundaryJac[nSpace*(nSpace-1)],
918  metricTensor[(nSpace-1)*(nSpace-1)],
919  metricTensorDetSqrt,
920  dS,
921  u_test_dS[nDOF_test_element],
922  u_grad_trial_trace[nDOF_trial_element*nSpace],
923  u_grad_test_dS[nDOF_trial_element*nSpace],
924  normal[nSpace],x_ext,y_ext,z_ext,xt_ext,yt_ext,zt_ext,integralScaling,
925  //
926  G[nSpace*nSpace],G_dd_G,tr_G;
927  //
928  //calculate the solution and gradients at quadrature points
929  //
930  //compute information about mapping from reference element to physical element
931  ck.calculateMapping_elementBoundary(eN,
932  ebN_local,
933  kb,
934  ebN_local_kb,
935  mesh_dof.data(),
936  mesh_l2g.data(),
937  mesh_trial_trace_ref.data(),
938  mesh_grad_trial_trace_ref.data(),
939  boundaryJac_ref.data(),
940  jac_ext,
941  jacDet_ext,
942  jacInv_ext,
943  boundaryJac,
944  metricTensor,
945  metricTensorDetSqrt,
946  normal_ref.data(),
947  normal,
948  x_ext,y_ext,z_ext);
949  /* ck.calculateMappingVelocity_elementBoundary(eN, */
950  /* ebN_local, */
951  /* kb, */
952  /* ebN_local_kb, */
953  /* mesh_velocity_dof, */
954  /* mesh_l2g, */
955  /* mesh_trial_trace_ref, */
956  /* xt_ext,yt_ext,zt_ext, */
957  /* normal, */
958  /* boundaryJac, */
959  /* metricTensor, */
960  /* integralScaling); */
961  /*dS = ((1.0-MOVING_DOMAIN)*metricTensorDetSqrt + MOVING_DOMAIN*integralScaling)*dS_ref.data()[kb];*/
962  dS = metricTensorDetSqrt*dS_ref.data()[kb];
963  //get the metric tensor
964  //cek todo use symmetry
965  ck.calculateG(jacInv_ext,G,G_dd_G,tr_G);
966  //compute shape and solution information
967  //shape
968  ck.gradTrialFromRef(&u_grad_trial_trace_ref.data()[ebN_local_kb_nSpace*nDOF_trial_element],jacInv_ext,u_grad_trial_trace);
969  //solution and gradients
970  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);
971  ck.gradFromDOF(u_dof.data(),&u_l2g.data()[eN_nDOF_trial_element],u_grad_trial_trace,grad_u_ext);
972  //precalculate test function products with integration weights
973  for (int j=0;j<nDOF_trial_element;j++)
974  {
975  u_test_dS[j] = u_test_trace_ref.data()[ebN_local_kb*nDOF_test_element+j]*dS;
976  for (int I=0;I<nSpace;I++)
977  u_grad_test_dS[j*nSpace+I] = u_grad_trial_trace[j*nSpace+I]*dS;//cek hack, using trial
978  }
979  //
980  //load the boundary values
981  //
982  bc_u_ext = isDOFBoundary_u.data()[ebNE_kb]*ebqe_bc_u_ext.data()[ebNE_kb]+(1-isDOFBoundary_u.data()[ebNE_kb])*u_ext;
983  //
984  //
985  //calculate the pde coefficients using the solution and the boundary values for the solution
986  //
987  a_ext = &ebqe_a.data()[ebNE_kb*sd_rowptr[nSpace]];
988  for (int I=0;I<nSpace;I++)
989  {
990  f_ext[I] = ebqe_v.data()[ebNE_kb*nSpace+I]*u_ext;
991  df_ext[I] = ebqe_v.data()[ebNE_kb*nSpace+I];
992  bc_f_ext[I] = ebqe_v.data()[ebNE_kb*nSpace+I]*bc_u_ext;
993  bc_df_ext[I] = ebqe_v.data()[ebNE_kb*nSpace+I];
994  }
995  /* evaluateCoefficients(&ebqe_velocity_ext.data()[ebNE_kb_nSpace], */
996  /* u_ext, */
997  /* //VRANS */
998  /* porosity_ext, */
999  /* // */
1000  /* m_ext, */
1001  /* dm_ext, */
1002  /* f_ext, */
1003  /* df_ext); */
1004  /* evaluateCoefficients(&ebqe_velocity_ext.data()[ebNE_kb_nSpace], */
1005  /* bc_u_ext, */
1006  /* //VRANS */
1007  /* porosity_ext, */
1008  /* // */
1009  /* bc_m_ext, */
1010  /* bc_dm_ext, */
1011  /* bc_f_ext, */
1012  /* bc_df_ext); */
1013  //
1014  //moving mesh
1015  //
1016  /* double velocity_ext[nSpace]; */
1017  /* double mesh_velocity[3]; */
1018  /* mesh_velocity[0] = xt_ext; */
1019  /* mesh_velocity[1] = yt_ext; */
1020  /* mesh_velocity[2] = zt_ext; */
1021  /* for (int I=0;I<nSpace;I++) */
1022  /* velocity_ext[I] = ebqe_velocity_ext.data()[ebNE_kb_nSpace+0] - MOVING_DOMAIN*mesh_velocity[I]; */
1023  //
1024  //calculate the numerical fluxes
1025  //
1026  exteriorNumericalDiffusiveFlux(sd_rowptr.data(),
1027  sd_colind.data(),
1028  isDOFBoundary_u.data()[ebNE_kb],
1029  isDiffusiveFluxBoundary_u.data()[ebNE_kb],
1030  normal,
1031  a_ext,
1032  bc_u_ext,
1033  ebqe_bc_flux_u_ext.data()[ebNE_kb],
1034  a_ext,
1035  grad_u_ext,
1036  u_ext,
1037  ebqe_penalty_ext.data()[ebNE_kb],
1038  flux_diff_ext);
1039  exteriorNumericalAdvectiveFlux(isDOFBoundary_u.data()[ebNE_kb],
1040  isAdvectiveFluxBoundary_u.data()[ebNE_kb],
1041  normal,
1042  bc_u_ext,
1043  ebqe_bc_flux_u_ext.data()[ebNE_kb],
1044  u_ext,
1045  df_ext,
1046  flux_advect_ext);
1047  //ebqe_flux.data()[ebNE_kb] = flux_ext;
1048  //
1049  //update residuals
1050  //
1051  for (int i=0;i<nDOF_test_element;i++)
1052  {
1053  //int ebNE_kb_i = ebNE_kb*nDOF_test_element+i;
1054  elementResidual_u[i] += ck.ExteriorElementBoundaryFlux(flux_diff_ext+flux_advect_ext,u_test_dS[i])+
1055  ck.ExteriorElementBoundaryDiffusionAdjoint(isDOFBoundary_u.data()[ebNE_kb],
1056  isDiffusiveFluxBoundary_u.data()[ebNE_kb],
1057  eb_adjoint_sigma,
1058  u_ext,
1059  bc_u_ext,
1060  normal,
1061  sd_rowptr.data(),
1062  sd_colind.data(),
1063  a_ext,
1064  &u_grad_test_dS[i*nSpace]);
1065  }//i
1066  }//kb
1067  //
1068  //update the element and global residual storage
1069  //
1070  for (int i=0;i<nDOF_test_element;i++)
1071  {
1072  int eN_i = eN*nDOF_test_element+i;
1073  globalResidual.data()[offset_u+stride_u*u_l2g.data()[eN_i]] += elementResidual_u[i];
1074  }//i
1075  }//ebNE
1076  }
1077 
1078  inline void calculateElementJacobian(//element
1079  xt::pyarray<double>& mesh_trial_ref,
1080  xt::pyarray<double>& mesh_grad_trial_ref,
1081  xt::pyarray<double>& mesh_dof,
1082  xt::pyarray<int>& mesh_l2g,
1083  xt::pyarray<double>& x_ref,
1084  xt::pyarray<double>& dV_ref,
1085  xt::pyarray<double>& u_trial_ref,
1086  xt::pyarray<double>& u_grad_trial_ref,
1087  xt::pyarray<double>& u_test_ref,
1088  xt::pyarray<double>& u_grad_test_ref,
1089  xt::pyarray<double>& elementDiameter,
1090  xt::pyarray<double>& elementBoundaryDiameter,
1091  xt::pyarray<double>& nodeDiametersArray,
1092  xt::pyarray<double>& cfl,
1093  double Ct_sge,
1094  double sc_uref,
1095  double sc_alpha,
1096  double useMetrics,
1097  //element boundary
1098  xt::pyarray<double>& mesh_trial_trace_ref,
1099  xt::pyarray<double>& mesh_grad_trial_trace_ref,
1100  xt::pyarray<double>& dS_ref,
1101  xt::pyarray<double>& u_trial_trace_ref,
1102  xt::pyarray<double>& u_grad_trial_trace_ref,
1103  xt::pyarray<double>& u_test_trace_ref,
1104  xt::pyarray<double>& u_grad_test_trace_ref,
1105  xt::pyarray<double>& normal_ref,
1106  xt::pyarray<double>& boundaryJac_ref,
1107  //physics
1108  int nElements_global,
1109  int nElementBoundaries_owned,
1110  xt::pyarray<int>& u_l2g,
1111  xt::pyarray<double>& u_dof,
1112  xt::pyarray<int>& sd_rowptr,
1113  xt::pyarray<int>& sd_colind,
1114  xt::pyarray<double>& q_a,
1115  xt::pyarray<double>& q_v,
1116  xt::pyarray<double>& q_r,
1117  int lag_shockCapturing,
1118  double shockCapturingDiffusion,
1119  xt::pyarray<double>& q_numDiff_u,
1120  xt::pyarray<double>& q_numDiff_u_last,
1121  xt::pyarray<double>& elementJacobian_u_u,
1122  xt::pyarray<double>& element_u,
1123  int eN,
1124  const bool embeddedBoundary,
1125  const double embeddedBoundary_penalty,
1126  xt::pyarray<double>& embeddedBoundary_normal_q,
1127  xt::pyarray<double>& embeddedBoundary_u_q)
1128  {
1129  for (int i=0;i<nDOF_test_element;i++)
1130  for (int j=0;j<nDOF_trial_element;j++)
1131  {
1132  elementJacobian_u_u.data()[i*nDOF_trial_element+j]=0.0;
1133  }
1134  for (int k=0;k<nQuadraturePoints_element;k++)
1135  {
1136  gf_s.set_quad(k);
1137  int eN_k = eN*nQuadraturePoints_element+k; //index to a scalar at a quadrature point
1138  int eN_k_3d = eN_k*3;
1139  //declare local storage
1140  register double u=0.0,
1141  grad_u[nSpace],
1142  m=0.0,dm=0.0,
1143  h_phi=0.0,
1144  r_s=0.0,dr_s=0.0,
1145  f[nSpace],df[nSpace],
1146  f_s[nSpace]={0.,0.},df_s[nSpace]={0.,0.},
1147  ham_s=0.0,dham_s[nSpace]={0.,0.},
1148  m_t=0.0,dm_t=0.0,
1149  dpdeResidual_u_u[nDOF_trial_element],
1150  Lstar_u[nDOF_test_element],
1151  dsubgridError_u_u[nDOF_trial_element],
1152  tau=0.0,tau0=0.0,tau1=0.0,
1153  *a=NULL,
1154  dr=0.0,
1155  jac[nSpace*nSpace],
1156  jacDet,
1157  jacInv[nSpace*nSpace],
1158  u_grad_trial[nDOF_trial_element*nSpace],
1159  dV,
1160  u_test_dV[nDOF_test_element],
1161  u_grad_test_dV[nDOF_test_element*nSpace],
1162  x,y,z,
1163  G[nSpace*nSpace],G_dd_G,tr_G;
1164  //
1165  //calculate solution and gradients at quadrature points
1166  //
1167  ck.calculateMapping_element(eN,
1168  k,
1169  mesh_dof.data(),
1170  mesh_l2g.data(),
1171  mesh_trial_ref.data(),
1172  mesh_grad_trial_ref.data(),
1173  jac,
1174  jacDet,
1175  jacInv,
1176  x,y,z);
1177  ck.calculateH_element(eN,
1178  k,
1179  nodeDiametersArray.data(),
1180  mesh_l2g.data(),
1181  mesh_trial_ref.data(),
1182  h_phi);
1183  //get the physical integration weight
1184  dV = fabs(jacDet)*dV_ref.data()[k];
1185  //get metric tensor and friends
1186  ck.calculateG(jacInv,G,G_dd_G,tr_G);
1187  //get the trial function gradients
1188  ck.gradTrialFromRef(&u_grad_trial_ref.data()[k*nDOF_trial_element*nSpace],jacInv,u_grad_trial);
1189  //get the solution
1190  ck.valFromElementDOF(element_u.data(),&u_trial_ref.data()[k*nDOF_trial_element],u);
1191  //get the solution gradients
1192  ck.gradFromElementDOF(element_u.data(),u_grad_trial,grad_u);
1193  //precalculate test function products with integration weights
1194  for (int j=0;j<nDOF_trial_element;j++)
1195  {
1196  u_test_dV[j] = u_test_ref.data()[k*nDOF_trial_element+j]*dV;
1197  for (int I=0;I<nSpace;I++)
1198  {
1199  u_grad_test_dV[j*nSpace+I] = u_grad_trial[j*nSpace+I]*dV;//cek warning won't work for Petrov-Galerkin
1200  }
1201  }
1202  //
1203  //calculate pde coefficients and derivatives at quadrature points
1204  //
1205  //evaluateCoefficients()
1206  a = &q_a.data()[eN_k*sd_rowptr.data()[nSpace]];
1207  for (int I=0;I<nSpace;I++)
1208  df[I] = q_v.data()[eN_k*nSpace+I];
1209  dr = 0.0;
1210  const double H_s = gf_s.H(0.,0.);
1211  const double D_s = gf_s.D(0.,0.);
1212  if(embeddedBoundary)
1213  {
1214  double level_set_normal[nSpace];
1215  double sign=0.0;
1216  double norm_exact=0.0,norm_cut=0.0;
1217  for (int I=0;I<nSpace;I++)
1218  {
1219  sign += embeddedBoundary_normal_q.data()[eN_k_3d+I]*gf_s.get_normal()[I];
1220  level_set_normal[I] = gf_s.get_normal()[I];
1221  norm_cut += level_set_normal[I]*level_set_normal[I];
1222  norm_exact += embeddedBoundary_normal_q.data()[eN_k_3d+I]*embeddedBoundary_normal_q.data()[eN_k_3d+I];
1223  }
1224  assert(std::fabs(1.0-norm_cut) < 1.0e-8);
1225  assert(std::fabs(1.0-norm_exact) < 1.0e-8);
1226  if (sign < 0.0)
1227  for (int I=0;I<nSpace;I++)
1228  level_set_normal[I]*=-1.0;
1229  updateEmbeddedBoundaryTerms(embeddedBoundary_penalty/h_phi,//penalty,
1230  dV,
1231  level_set_normal,
1232  embeddedBoundary_u_q.data()[eN_k],
1233  u,
1234  grad_u,
1235  a[0],//assume scalar diffusion for now
1236  r_s,
1237  dr_s,
1238  ham_s,
1239  dham_s,
1240  f_s,
1241  df_s);
1242  }
1243  //
1244  //calculate subgrid error contribution to the Jacobian (strong residual, adjoint, jacobian of strong residual)
1245  //
1246  //calculate the adjoint times the test functions
1247  for (int i=0;i<nDOF_test_element;i++)
1248  {
1249  // int eN_k_i_nSpace = (eN_k*nDOF_trial_element+i)*nSpace;
1250  // Lstar_u[i]=ck.Advection_adjoint(df,&u_grad_test_dV.data()[eN_k_i_nSpace]);
1251  register int i_nSpace = i*nSpace;
1252  Lstar_u[i]=ck.Advection_adjoint(df,&u_grad_test_dV[i_nSpace]);
1253  }
1254  //calculate the Jacobian of strong residual
1255  for (int j=0;j<nDOF_trial_element;j++)
1256  {
1257  //int eN_k_j=eN_k*nDOF_trial_element+j;
1258  //int eN_k_j_nSpace = eN_k_j*nSpace;
1259  int j_nSpace = j*nSpace;
1260  dpdeResidual_u_u[j]= ck.MassJacobian_strong(dm_t,u_trial_ref.data()[k*nDOF_trial_element+j]) +
1261  ck.AdvectionJacobian_strong(df,&u_grad_trial[j_nSpace]);
1262  }
1263  //tau and tau*Res
1264  calculateSubgridError_tau(elementDiameter.data()[eN],
1265  dm_t,
1266  df,
1267  cfl.data()[eN_k],
1268  tau0);
1269 
1270  calculateSubgridError_tau(Ct_sge,
1271  G,
1272  dm_t,
1273  df,
1274  tau1,
1275  cfl.data()[eN_k]);
1276  tau = useMetrics*tau1+(1.0-useMetrics)*tau0;
1277 
1278  for(int j=0;j<nDOF_trial_element;j++)
1279  dsubgridError_u_u[j] = -tau*dpdeResidual_u_u[j];
1280 
1281  for(int i=0;i<nDOF_test_element;i++)
1282  {
1283  int i_nSpace=i*nSpace;
1284  for(int j=0;j<nDOF_trial_element;j++)
1285  {
1286  int j_nSpace = j*nSpace;
1287  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]) +
1288  ck.SimpleDiffusionJacobian_weak(sd_rowptr.data(),sd_colind.data(),a,&u_grad_trial[j_nSpace],&u_grad_test_dV[i_nSpace]) +
1289  ck.ReactionJacobian_weak(dr,u_trial_ref.data()[k*nDOF_trial_element+j],u_test_dV[i]) +
1290  ck.SubgridErrorJacobian(dsubgridError_u_u[j],Lstar_u[i]) +
1291  ck.NumericalDiffusionJacobian(q_numDiff_u_last.data()[eN_k],&u_grad_trial[j_nSpace],&u_grad_test_dV[i_nSpace]));
1292  if (embeddedBoundary)
1293  {
1294  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]) +
1295  ck.ReactionJacobian_weak(dr_s,u_trial_ref.data()[k*nDOF_trial_element+j], u_test_dV[i]) +
1296  ck.HamiltonianJacobian_weak(dham_s,&u_grad_trial[j_nSpace],u_test_dV[i]));
1297  }
1298  }//j
1299  }//i
1300  }//k
1301  }
1302 
1303  void calculateJacobian(arguments_dict& args)
1304  {
1305  gf.useExact=true;
1306  gf_s.useExact=true;
1307  xt::pyarray<double>& mesh_trial_ref = args.array<double>("mesh_trial_ref");
1308  xt::pyarray<double>& mesh_grad_trial_ref = args.array<double>("mesh_grad_trial_ref");
1309  xt::pyarray<double>& mesh_dof = args.array<double>("mesh_dof");
1310  xt::pyarray<int>& mesh_l2g = args.array<int>("mesh_l2g");
1311  xt::pyarray<double>& dV_ref = args.array<double>("dV_ref");
1312  xt::pyarray<double>& u_trial_ref = args.array<double>("u_trial_ref");
1313  xt::pyarray<double>& u_grad_trial_ref = args.array<double>("u_grad_trial_ref");
1314  xt::pyarray<double>& u_test_ref = args.array<double>("u_test_ref");
1315  xt::pyarray<double>& u_grad_test_ref = args.array<double>("u_grad_test_ref");
1316  xt::pyarray<double>& elementDiameter = args.array<double>("elementDiameter");
1317  xt::pyarray<double>& cfl = args.array<double>("cfl");
1318  double Ct_sge = args.scalar<double>("Ct_sge");
1319  double sc_uref = args.scalar<double>("sc_uref");
1320  double sc_alpha = args.scalar<double>("sc_alpha");
1321  double useMetrics = args.scalar<double>("useMetrics");
1322  xt::pyarray<double>& mesh_trial_trace_ref = args.array<double>("mesh_trial_trace_ref");
1323  xt::pyarray<double>& mesh_grad_trial_trace_ref = args.array<double>("mesh_grad_trial_trace_ref");
1324  xt::pyarray<double>& dS_ref = args.array<double>("dS_ref");
1325  xt::pyarray<double>& u_trial_trace_ref = args.array<double>("u_trial_trace_ref");
1326  xt::pyarray<double>& u_grad_trial_trace_ref = args.array<double>("u_grad_trial_trace_ref");
1327  xt::pyarray<double>& u_test_trace_ref = args.array<double>("u_test_trace_ref");
1328  xt::pyarray<double>& u_grad_test_trace_ref = args.array<double>("u_grad_test_trace_ref");
1329  xt::pyarray<double>& normal_ref = args.array<double>("normal_ref");
1330  xt::pyarray<double>& boundaryJac_ref = args.array<double>("boundaryJac_ref");
1331  int nElements_global = args.scalar<int>("nElements_global");
1332  xt::pyarray<int>& u_l2g = args.array<int>("u_l2g");
1333  xt::pyarray<double>& u_dof = args.array<double>("u_dof");
1334  xt::pyarray<int>& sd_rowptr = args.array<int>("sd_rowptr");
1335  xt::pyarray<int>& sd_colind = args.array<int>("sd_colind");
1336  xt::pyarray<double>& q_a = args.array<double>("q_a");
1337  xt::pyarray<double>& q_v = args.array<double>("q_v");
1338  xt::pyarray<double>& q_r = args.array<double>("q_r");
1339  int lag_shockCapturing = args.scalar<int>("lag_shockCapturing");
1340  double shockCapturingDiffusion = args.scalar<double>("shockCapturingDiffusion");
1341  xt::pyarray<double>& q_numDiff_u = args.array<double>("q_numDiff_u");
1342  xt::pyarray<double>& q_numDiff_u_last = args.array<double>("q_numDiff_u_last");
1343  xt::pyarray<int>& csrRowIndeces_u_u = args.array<int>("csrRowIndeces_u_u");
1344  xt::pyarray<int>& csrColumnOffsets_u_u = args.array<int>("csrColumnOffsets_u_u");
1345  xt::pyarray<double>& globalJacobian = args.array<double>("globalJacobian");
1346  int nExteriorElementBoundaries_global = args.scalar<int>("nExteriorElementBoundaries_global");
1347  xt::pyarray<int>& exteriorElementBoundariesArray = args.array<int>("exteriorElementBoundariesArray");
1348  xt::pyarray<int>& elementBoundaryElementsArray = args.array<int>("elementBoundaryElementsArray");
1349  xt::pyarray<int>& elementBoundaryLocalElementBoundariesArray = args.array<int>("elementBoundaryLocalElementBoundariesArray");
1350  xt::pyarray<double>& ebqe_a = args.array<double>("ebqe_a");
1351  xt::pyarray<double>& ebqe_v = args.array<double>("ebqe_v");
1352  xt::pyarray<int>& isDOFBoundary_u = args.array<int>("isDOFBoundary_u");
1353  xt::pyarray<double>& ebqe_bc_u_ext = args.array<double>("ebqe_bc_u_ext");
1354  xt::pyarray<int>& isDiffusiveFluxBoundary_u = args.array<int>("isDiffusiveFluxBoundary_u");
1355  xt::pyarray<int>& isAdvectiveFluxBoundary_u = args.array<int>("isAdvectiveFluxBoundary_u");
1356  xt::pyarray<double>& ebqe_bc_flux_u_ext = args.array<double>("ebqe_bc_flux_u_ext");
1357  xt::pyarray<double>& ebqe_bc_advectiveFlux_u_ext = args.array<double>("ebqe_bc_advectiveFlux_u_ext");
1358  xt::pyarray<int>& csrColumnOffsets_eb_u_u = args.array<int>("csrColumnOffsets_eb_u_u");
1359  xt::pyarray<double>& ebqe_penalty_ext = args.array<double>("ebqe_penalty_ext");
1360  const bool embeddedBoundary = args.scalar<int>("embeddedBoundary");
1361  const double embeddedBoundary_penalty = args.scalar<double>("embeddedBoundary_penalty");
1362  const double embeddedBoundary_ghost_penalty = args.scalar<double>("embedded_ghost_penalty");
1363  xt::pyarray<double>& embeddedBoundary_sdf_nodes = args.array<double>("embeddedBoundary_sdf_nodes");
1364  xt::pyarray<double>& embeddedBoundary_sdf_q = args.array<double>("embeddedBoundary_sdf_q");
1365  xt::pyarray<double>& embeddedBoundary_normal_q = args.array<double>("embeddedBoundary_normal_q");
1366  xt::pyarray<double>& embeddedBoundary_u_q = args.array<double>("embeddedBoundary_u_q");
1367  xt::pyarray<double>& isActiveDOF = args.array<double>("isActiveDOF");
1368  const double eb_adjoint_sigma = args.scalar<double>("eb_adjoint_sigma");
1369  xt::pyarray<double>& x_ref = args.array<double>("x_ref");
1370  xt::pyarray<int>& elementBoundariesArray = args.array<int>("elementBoundariesArray");
1371  const int nElementBoundaries_owned = args.scalar<int>("nElementBoundaries_owned");
1372  xt::pyarray<double>& elementBoundaryDiameter = args.array<double>("elementBoundaryDiameter");
1373  xt::pyarray<double>& nodeDiametersArray = args.array<double>("nodeDiametersArray");
1374  //
1375  //loop over elements to compute volume integrals and load them into the element Jacobians and global Jacobian
1376  //
1377  for(int eN=0;eN<nElements_global;eN++)
1378  {
1379  //register double elementJacobian_u_u.data()[nDOF_test_element*nDOF_trial_element],element_u[nDOF_trial_element];
1380  auto elementJacobian_u_u = xt::pyarray<double>::from_shape({nDOF_test_element*nDOF_trial_element});
1381  auto element_u = xt::pyarray<double>::from_shape({nDOF_trial_element});
1382  for (int j=0;j<nDOF_trial_element;j++)
1383  {
1384  register int eN_j = eN*nDOF_trial_element+j;
1385  element_u.data()[j] = u_dof.data()[u_l2g.data()[eN_j]];
1386  }
1387  double element_phi_s[nDOF_mesh_trial_element];
1388  for (int j=0;j<nDOF_mesh_trial_element;j++)
1389  {
1390  register int eN_j = eN*nDOF_mesh_trial_element+j;
1391  element_phi_s[j] = embeddedBoundary_sdf_nodes.data()[u_l2g.data()[eN_j]];
1392  }
1393  double element_nodes[nDOF_mesh_trial_element*3];
1394  for (int i=0;i<nDOF_mesh_trial_element;i++)
1395  {
1396  register int eN_i=eN*nDOF_mesh_trial_element+i;
1397  for(int I=0;I<3;I++)
1398  element_nodes[i*3 + I] = mesh_dof.data()[mesh_l2g.data()[eN_i]*3 + I];
1399  }//i
1400  int icase_s = gf_s.calculate(element_phi_s, element_nodes, x_ref.data(), false);
1401  //int icase = gf.calculate(element_phi, element_nodes, x_ref.data(), rho_1*nu_1, rho_0*nu_0,false,false);
1402  calculateElementJacobian(mesh_trial_ref,
1403  mesh_grad_trial_ref,
1404  mesh_dof,
1405  mesh_l2g,
1406  x_ref,
1407  dV_ref,
1408  u_trial_ref,
1409  u_grad_trial_ref,
1410  u_test_ref,
1411  u_grad_test_ref,
1412  elementDiameter,
1413  elementBoundaryDiameter,
1414  nodeDiametersArray,
1415  cfl,
1416  Ct_sge,
1417  sc_uref,
1418  sc_alpha,
1419  useMetrics,
1420  mesh_trial_trace_ref,
1421  mesh_grad_trial_trace_ref,
1422  dS_ref,
1423  u_trial_trace_ref,
1424  u_grad_trial_trace_ref,
1425  u_test_trace_ref,
1426  u_grad_test_trace_ref,
1427  normal_ref,
1428  boundaryJac_ref,
1429  nElements_global,
1430  nElementBoundaries_owned,
1431  u_l2g,
1432  u_dof,
1433  sd_rowptr,
1434  sd_colind,
1435  q_a,
1436  q_v,
1437  q_r,
1438  lag_shockCapturing,
1439  shockCapturingDiffusion,
1440  q_numDiff_u,
1441  q_numDiff_u_last,
1442  elementJacobian_u_u,
1443  element_u,
1444  eN,
1445  embeddedBoundary,
1446  embeddedBoundary_penalty,
1447  embeddedBoundary_normal_q,
1448  embeddedBoundary_u_q);
1449  //
1450  //load into element Jacobian into global Jacobian
1451  //
1452  for (int i=0;i<nDOF_test_element;i++)
1453  {
1454  int eN_i = eN*nDOF_test_element+i;
1455  for (int j=0;j<nDOF_trial_element;j++)
1456  {
1457  int eN_i_j = eN_i*nDOF_trial_element+j;
1458  globalJacobian.data()[csrRowIndeces_u_u.data()[eN_i] + csrColumnOffsets_u_u.data()[eN_i_j]] += elementJacobian_u_u.data()[i*nDOF_trial_element+j];
1459  }//j
1460  }//i
1461  }//elements
1462  for (std::set<int>::iterator it=cutfem_boundaries.begin(); it!=cutfem_boundaries.end(); ++it)
1463  {
1464  std::map<int,double> Dw_Dn_jump;
1465  std::map<std::pair<int, int>, int> u_u_nz;
1466  register double gamma_cutfem=embeddedBoundary_ghost_penalty,h_cutfem=elementBoundaryDiameter.data()[*it];
1467  int eN_nDOF_trial_element = elementBoundaryElementsArray.data()[(*it)*2+0]*nDOF_trial_element;
1468  for (int kb=0;kb<nQuadraturePoints_elementBoundary;kb++)
1469  {
1470  register double Dp_Dn_jump=0.0, Du_Dn_jump=0.0, Dv_Dn_jump=0.0,dS;
1471  for (int eN_side=0;eN_side < 2; eN_side++)
1472  {
1473  register int ebN = *it,
1474  eN = elementBoundaryElementsArray.data()[ebN*2+eN_side];
1475  for (int i=0;i<nDOF_test_element;i++)
1476  Dw_Dn_jump[u_l2g.data()[eN*nDOF_test_element+i]] = 0.0;
1477  }
1478  for (int eN_side=0;eN_side < 2; eN_side++)
1479  {
1480  register int ebN = *it,
1481  eN = elementBoundaryElementsArray.data()[ebN*2+eN_side],
1482  ebN_local = elementBoundaryLocalElementBoundariesArray.data()[ebN*2+eN_side],
1483  eN_nDOF_trial_element = eN*nDOF_trial_element,
1484  ebN_local_kb = ebN_local*nQuadraturePoints_elementBoundary+kb,
1485  ebN_local_kb_nSpace = ebN_local_kb*nSpace;
1486  register double
1487  u_int=0.0,
1488  grad_u_int[nSpace]={0.,0.},
1489  jac_int[nSpace*nSpace],
1490  jacDet_int,
1491  jacInv_int[nSpace*nSpace],
1492  boundaryJac[nSpace*(nSpace-1)],
1493  metricTensor[(nSpace-1)*(nSpace-1)],
1494  metricTensorDetSqrt,
1495  u_test_dS[nDOF_test_element],
1496  u_grad_trial_trace[nDOF_trial_element*nSpace],
1497  u_grad_test_dS[nDOF_trial_element*nSpace],
1498  normal[2],x_int,y_int,z_int,xt_int,yt_int,zt_int,integralScaling,
1499  G[nSpace*nSpace],G_dd_G,tr_G,h_phi,h_penalty,penalty;
1500  //compute information about mapping from reference element to physical element
1501  ck.calculateMapping_elementBoundary(eN,
1502  ebN_local,
1503  kb,
1504  ebN_local_kb,
1505  mesh_dof.data(),
1506  mesh_l2g.data(),
1507  mesh_trial_trace_ref.data(),
1508  mesh_grad_trial_trace_ref.data(),
1509  boundaryJac_ref.data(),
1510  jac_int,
1511  jacDet_int,
1512  jacInv_int,
1513  boundaryJac,
1514  metricTensor,
1515  metricTensorDetSqrt,
1516  normal_ref.data(),
1517  normal,
1518  x_int,y_int,z_int);
1519  dS = metricTensorDetSqrt*dS_ref.data()[kb];
1520  //compute shape and solution information
1521  //shape
1522  ck.gradTrialFromRef(&u_grad_trial_trace_ref.data()[ebN_local_kb_nSpace*nDOF_trial_element],jacInv_int,u_grad_trial_trace);
1523  for (int i=0;i<nDOF_test_element;i++)
1524  {
1525  int eN_i = eN*nDOF_test_element + i;
1526  for (int I=0;I<nSpace;I++)
1527  Dw_Dn_jump[u_l2g.data()[eN_i]] += u_grad_trial_trace[i*nSpace+I]*normal[I];
1528  }
1529  }//eN_side
1530  for (int eN_side=0;eN_side < 2; eN_side++)
1531  {
1532  register int ebN = *it,
1533  eN = elementBoundaryElementsArray.data()[ebN*2+eN_side];
1534  for (int i=0;i<nDOF_test_element;i++)
1535  {
1536  register int eN_i = eN*nDOF_test_element+i;
1537  for (int eN_side2=0;eN_side2 < 2; eN_side2++)
1538  {
1539  register int eN2 = elementBoundaryElementsArray.data()[ebN*2+eN_side2];
1540  for (int j=0;j<nDOF_test_element;j++)
1541  {
1542  int eN_i_j = eN_i*nDOF_test_element + j;
1543  int eN2_j = eN2*nDOF_test_element + j;
1544  register int ebN_i_j = ebN*4*nDOF_test_X_trial_element +
1545  eN_side*2*nDOF_test_X_trial_element +
1546  eN_side2*nDOF_test_X_trial_element +
1547  i*nDOF_trial_element +
1548  j;
1549  std::pair<int,int> ij = std::make_pair(u_l2g.data()[eN_i], u_l2g.data()[eN2_j]);
1550  if (u_u_nz.count(ij))
1551  {
1552  assert(u_u_nz[ij] == csrRowIndeces_u_u.data()[eN_i] + csrColumnOffsets_eb_u_u.data()[ebN_i_j]);
1553  }
1554  else
1555  u_u_nz[ij] = csrRowIndeces_u_u.data()[eN_i] + csrColumnOffsets_eb_u_u.data()[ebN_i_j];
1556  }
1557  }
1558  }
1559  }
1560  for (std::map<int,double>::iterator wi_it=Dw_Dn_jump.begin(); wi_it!=Dw_Dn_jump.end(); ++wi_it)
1561  for (std::map<int,double>::iterator wj_it=Dw_Dn_jump.begin(); wj_it!=Dw_Dn_jump.end(); ++wj_it)
1562  {
1563  int i_global = wi_it->first,
1564  j_global = wj_it->first;
1565  double Dw_Dn_jump_i = wi_it->second,
1566  Dw_Dn_jump_j = wj_it->second;
1567  std::pair<int,int> ij = std::make_pair(i_global, j_global);
1568  globalJacobian.data()[u_u_nz.at(ij)] += gamma_cutfem*h_cutfem*Dw_Dn_jump_j*Dw_Dn_jump_i*dS;
1569  }//i,j
1570  }//kb
1571  }//cutfem element boundaries
1572  //
1573  //loop over exterior element boundaries to compute the surface integrals and load them into the global Jacobian
1574  //
1575  for (int ebNE = 0; ebNE < nExteriorElementBoundaries_global; ebNE++)
1576  {
1577  register int ebN = exteriorElementBoundariesArray.data()[ebNE];
1578  register int eN = elementBoundaryElementsArray.data()[ebN*2+0],
1579  ebN_local = elementBoundaryLocalElementBoundariesArray.data()[ebN*2+0],
1580  eN_nDOF_trial_element = eN*nDOF_trial_element;
1581  for (int kb=0;kb<nQuadraturePoints_elementBoundary;kb++)
1582  {
1583  register int ebNE_kb = ebNE*nQuadraturePoints_elementBoundary+kb,
1584  ebNE_kb_nSpace = ebNE_kb*nSpace,
1585  ebN_local_kb = ebN_local*nQuadraturePoints_elementBoundary+kb,
1586  ebN_local_kb_nSpace = ebN_local_kb*nSpace;
1587 
1588  register double u_ext=0.0,
1589  grad_u_ext[nSpace],
1590  m_ext=0.0,
1591  dm_ext=0.0,
1592  *a_ext,
1593  f_ext[nSpace],
1594  df_ext[nSpace],
1595  r_ext=0.0,
1596  dflux_u_u_ext=0.0,
1597  bc_u_ext=0.0,
1598  //bc_grad_u_ext[nSpace],
1599  bc_m_ext=0.0,
1600  bc_dm_ext=0.0,
1601  bc_f_ext[nSpace],
1602  bc_df_ext[nSpace],
1603  fluxJacobian_u_u[nDOF_trial_element],
1604  jac_ext[nSpace*nSpace],
1605  jacDet_ext,
1606  jacInv_ext[nSpace*nSpace],
1607  boundaryJac[nSpace*(nSpace-1)],
1608  metricTensor[(nSpace-1)*(nSpace-1)],
1609  metricTensorDetSqrt,
1610  dS,
1611  u_test_dS[nDOF_test_element],
1612  u_grad_trial_trace[nDOF_trial_element*nSpace],
1613  u_grad_test_dS[nDOF_trial_element*nSpace],
1614  normal[nSpace],x_ext,y_ext,z_ext,xt_ext,yt_ext,zt_ext,integralScaling,
1615  //
1616  G[nSpace*nSpace],G_dd_G,tr_G;
1617  ck.calculateMapping_elementBoundary(eN,
1618  ebN_local,
1619  kb,
1620  ebN_local_kb,
1621  mesh_dof.data(),
1622  mesh_l2g.data(),
1623  mesh_trial_trace_ref.data(),
1624  mesh_grad_trial_trace_ref.data(),
1625  boundaryJac_ref.data(),
1626  jac_ext,
1627  jacDet_ext,
1628  jacInv_ext,
1629  boundaryJac,
1630  metricTensor,
1631  metricTensorDetSqrt,
1632  normal_ref.data(),
1633  normal,
1634  x_ext,y_ext,z_ext);
1635  dS = metricTensorDetSqrt*dS_ref.data()[kb];
1636  ck.calculateG(jacInv_ext,G,G_dd_G,tr_G);
1637  //compute shape and solution information
1638  //shape
1639  ck.gradTrialFromRef(&u_grad_trial_trace_ref.data()[ebN_local_kb_nSpace*nDOF_trial_element],jacInv_ext,u_grad_trial_trace);
1640  //solution and gradients
1641  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);
1642  ck.gradFromDOF(u_dof.data(),&u_l2g.data()[eN_nDOF_trial_element],u_grad_trial_trace,grad_u_ext);
1643  //precalculate test function products with integration weights
1644  for (int j=0;j<nDOF_trial_element;j++)
1645  {
1646  u_test_dS[j] = u_test_trace_ref.data()[ebN_local_kb*nDOF_test_element+j]*dS;
1647  for (int I=0;I<nSpace;I++)
1648  u_grad_test_dS[j*nSpace+I] = u_grad_trial_trace[j*nSpace+I]*dS;//cek hack, using trial
1649  }
1650  //
1651  //load the boundary values
1652  //
1653  bc_u_ext = isDOFBoundary_u.data()[ebNE_kb]*ebqe_bc_u_ext.data()[ebNE_kb]+(1-isDOFBoundary_u.data()[ebNE_kb])*u_ext;
1654  a_ext = &ebqe_a.data()[ebNE_kb*sd_rowptr.data()[nSpace]];
1655  for (int I=0;I<nSpace;I++)
1656  {
1657  df_ext[I] = ebqe_v.data()[ebNE_kb*nSpace+I];
1658  bc_df_ext[I] = ebqe_v.data()[ebNE_kb*nSpace+I];
1659  }
1660  //
1661  //calculate the numerical fluxes
1662  //
1663  exteriorNumericalAdvectiveFluxDerivative(isDOFBoundary_u.data()[ebNE_kb],
1664  isAdvectiveFluxBoundary_u.data()[ebNE_kb],
1665  normal,
1666  df_ext,
1667  dflux_u_u_ext);
1668  //
1669  //calculate the flux jacobian
1670  //
1671  for (int j=0;j<nDOF_trial_element;j++)
1672  {
1673  //register int ebNE_kb_j = ebNE_kb*nDOF_trial_element+j;
1674  register int j_nSpace = j*nSpace,ebN_local_kb_j=ebN_local_kb*nDOF_trial_element+j;
1675  fluxJacobian_u_u[j]=ExteriorNumericalDiffusiveFluxJacobian(sd_rowptr.data(),
1676  sd_colind.data(),
1677  isDOFBoundary_u.data()[ebNE_kb],
1678  isDiffusiveFluxBoundary_u.data()[ebNE_kb],
1679  normal,
1680  a_ext,
1681  u_trial_trace_ref.data()[ebN_local_kb_j],
1682  &u_grad_trial_trace[j_nSpace],
1683  ebqe_penalty_ext.data()[ebNE_kb]) +
1684  ck.ExteriorNumericalAdvectiveFluxJacobian(dflux_u_u_ext,u_trial_trace_ref.data()[ebN_local_kb_j]);
1685  }//j
1686  //
1687  //update the global Jacobian from the flux Jacobian
1688  //
1689  for (int i=0;i<nDOF_test_element;i++)
1690  {
1691  register int eN_i = eN*nDOF_test_element+i;
1692  register int i_nSpace = i*nSpace;
1693  for (int j=0;j<nDOF_trial_element;j++)
1694  {
1695  register int ebN_i_j = ebN*4*nDOF_test_X_trial_element + i*nDOF_trial_element + j;
1696  register int ebN_local_kb_j=ebN_local_kb*nDOF_trial_element+j;
1697 
1698  globalJacobian.data()[csrRowIndeces_u_u.data()[eN_i] + csrColumnOffsets_eb_u_u.data()[ebN_i_j]] += fluxJacobian_u_u[j]*u_test_dS[i]+
1699  ck.ExteriorElementBoundaryDiffusionAdjointJacobian(isDOFBoundary_u.data()[ebNE_kb],
1700  isDiffusiveFluxBoundary_u.data()[ebNE_kb],
1701  eb_adjoint_sigma,
1702  u_trial_trace_ref.data()[ebN_local_kb_j],
1703  normal,
1704  sd_rowptr.data(),
1705  sd_colind.data(),
1706  a_ext,
1707  &u_grad_test_dS[i_nSpace]);
1708  }//j
1709  }//i
1710  }//kb
1711  }//ebNE
1712  }//computeJacobian
1713  };//cADR
1714 
1715  inline cADR_base* newADR(int nSpaceIn,
1716  int nQuadraturePoints_elementIn,
1717  int nDOF_mesh_trial_elementIn,
1718  int nDOF_trial_elementIn,
1719  int nDOF_test_elementIn,
1720  int nQuadraturePoints_elementBoundaryIn,
1721  int CompKernelFlag)
1722  {
1723  if (nSpaceIn == 2)
1724  return proteus::chooseAndAllocateDiscretization2D<cADR_base,cADR,CompKernel>(nSpaceIn,
1725  nQuadraturePoints_elementIn,
1726  nDOF_mesh_trial_elementIn,
1727  nDOF_trial_elementIn,
1728  nDOF_test_elementIn,
1729  nQuadraturePoints_elementBoundaryIn,
1730  CompKernelFlag);
1731  else
1732  return proteus::chooseAndAllocateDiscretization<cADR_base,cADR,CompKernel>(nSpaceIn,
1733  nQuadraturePoints_elementIn,
1734  nDOF_mesh_trial_elementIn,
1735  nDOF_trial_elementIn,
1736  nDOF_test_elementIn,
1737  nQuadraturePoints_elementBoundaryIn,
1738  CompKernelFlag);
1739  }
1740 }//proteus
1741 #endif
sign
#define sign(x, y)
Definition: jf.h:44
proteus::cADR::calculateSubgridError_tau
void calculateSubgridError_tau(const double &elementDiameter, const double &dmt, const double dH[nSpace], double &cfl, double &tau)
Definition: ADR.h:133
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:153
proteus::cADR::ifem_boundaries
std::set< int > ifem_boundaries
Definition: ADR.h:40
proteus::cADR_base::~cADR_base
virtual ~cADR_base()
Definition: ADR.h:24
proteus::cADR::ck
CompKernelType ck
Definition: ADR.h:44
proteus::cADR_base::calculateJacobian
virtual void calculateJacobian(arguments_dict &args)=0
n
Int n
Definition: Headers.h:28
df
double df(double C, double b, double a, int q, int r)
Definition: analyticalSolutions.c:2209
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)
Definition: ADR.h:271
ModelFactory.h
CompKernel.h
proteus::arguments_dict::scalar
T & scalar(const std::string &key)
proteus::arguments_dict::array
xt::pyarray< T > & array(const std::string &key)
num
Int num
Definition: Headers.h:32
proteus::cADR::nDOF_test_X_trial_element
const int nDOF_test_X_trial_element
Definition: ADR.h:43
proteus::cADR::calculateResidual
void calculateResidual(arguments_dict &args)
Definition: ADR.h:579
v
Double v
Definition: Headers.h:95
proteus::cADR::cutfem_boundary_elements
std::set< int > cutfem_boundary_elements
Definition: ADR.h:41
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:56
equivalent_polynomials.h
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:188
z
Double * z
Definition: Headers.h:49
u
Double u
Definition: Headers.h:89
proteus::cADR::ifem_boundary_elements
std::set< int > ifem_boundary_elements
Definition: ADR.h:40
proteus::cADR::calculateElementJacobian
void calculateElementJacobian(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)
Definition: ADR.h:1078
proteus::cADR_base::calculateResidual
virtual void calculateResidual(arguments_dict &args)=0
proteus::cADR::calculateJacobian
void calculateJacobian(arguments_dict &args)
Definition: ADR.h:1303
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:1715
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:102
proteus::cADR::gf
GeneralizedFunctions< nSpace, 3, nQuadraturePoints_element, nQuadraturePoints_elementBoundary > gf
Definition: ADR.h:45
equivalent_polynomials::GeneralizedFunctions_mix
Definition: equivalent_polynomials.h:767
proteus::cADR
Definition: ADR.h:38
proteus
Definition: ADR.h:17
proteus::cADR::gf_s
GeneralizedFunctions< nSpace, 3, nQuadraturePoints_element, nQuadraturePoints_elementBoundary > gf_s
Definition: ADR.h:46
proteus::f
double f(const double &g, const double &h, const double &hZ)
Definition: SW2DCV.h:58
r
Double r
Definition: Headers.h:83
proteus::arguments_dict
Definition: ArgumentsDict.h:70
proteus::cADR_base
Definition: ADR.h:22
proteus::cADR::calculateElementResidual
void calculateElementResidual(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, bool &element_active, std::valarray< bool > &elementIsActive)
Definition: ADR.h:307
proteus::cADR::cutfem_boundaries
std::set< int > cutfem_boundaries
Definition: ADR.h:41
proteus::cADR::cADR
cADR()
Definition: ADR.h:47
ArgumentsDict.h
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:234
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:168
proteus::cADR::elementIsActive
std::valarray< bool > elementIsActive
Definition: ADR.h:42
Generated on Fri Jul 1 2022 11:15:12 for proteus by doxygen 1.8.20