proteus  1.8.1
C/C++/Fortran libraries
RDLS.h
Go to the documentation of this file.
1 #ifndef RDLS_H
2 #define RDLS_H
3 #include <cmath>
4 #include <iostream>
5 #include <valarray>
6 #include "CompKernel.h"
7 #include "ModelFactory.h"
8 #include "equivalent_polynomials.h"
9 #include "ArgumentsDict.h"
10 #include "xtensor-python/pyarray.hpp"
11 
12 namespace py = pybind11;
13 
14 #define SINGLE_POTENTIAL 1
15 
16 namespace proteus
17 {
18  inline double heaviside(const double &z)
19  {
20  return (z>0 ? 1. : (z<0 ? 0. : 0.5));
21  }
22 
23  template<int nSpace, int nP, int nQ, int nEBQ>
24  using GeneralizedFunctions = equivalent_polynomials::GeneralizedFunctions_mix<nSpace, nP, nQ, nEBQ>;
25 
26 
27  class RDLS_base
28  {
29  public:
30  std::valarray<double> weighted_lumped_mass_matrix;
31  virtual ~RDLS_base(){}
32  virtual void calculateResidual(arguments_dict& args, bool useExact)=0;
33  virtual void calculateJacobian(arguments_dict& args, bool useExact)=0;
34  virtual void calculateResidual_ellipticRedist(arguments_dict& args, bool useExact)=0;
35  virtual void calculateJacobian_ellipticRedist(arguments_dict& args, bool useExact)=0;
36  virtual void normalReconstruction(arguments_dict& args)=0;
37  virtual std::tuple<double, double, double> calculateMetricsAtEOS(arguments_dict& args)=0;
38  };
39 
40  template<class CompKernelType,
41  int nSpace,
42  int nQuadraturePoints_element,
43  int nDOF_mesh_trial_element,
44  int nDOF_trial_element,
45  int nDOF_test_element,
46  int nQuadraturePoints_elementBoundary>
47  class RDLS : public RDLS_base
48  {
49  public:
50  const int nDOF_test_X_trial_element;
51  CompKernelType ck;
52  GeneralizedFunctions<nSpace,2,nQuadraturePoints_element,nQuadraturePoints_elementBoundary> gf,gfu;
53  RDLS():
54  nDOF_test_X_trial_element(nDOF_test_element*nDOF_trial_element),
55  ck()
56  {}
57 
58  inline
59  void evaluateCoefficients(const double& eps,
60  const double& u_levelSet,
61  const double& u,
62  const double grad_u[nSpace],
63  double& m,
64  double& dm,
65  double& H,
66  double dH[nSpace],
67  double& r)
68  {
69  int I;
70  double normGradU=0.0,Si=0.0;
71  m = u;
72  dm=1.0;
73  H = 0.0;
74  Si= gf.H(eps,u_levelSet) - gf.ImH(eps,u_levelSet);
75  /* if (u_levelSet > 0.0) */
76  /* Si=1.0; */
77  /* else if (u_levelSet < 0.0) */
78  /* Si = -1.0; */
79  /* else */
80  /* Si=0.0; */
81  r = -Si;
82  for (I=0; I < nSpace; I++)
83  {
84  normGradU += grad_u[I]*grad_u[I];
85  }
86  normGradU = sqrt(normGradU);
87  H = Si*normGradU;
88  for (I=0; I < nSpace; I++)
89  {
90  dH[I] = Si*grad_u[I]/(normGradU+1.0e-12);
91  }
92  }
93 
94  inline
95  void calculateSubgridError_tau(const double& elementDiameter,
96  const double& dmt,
97  const double dH[nSpace],
98  double& cfl,
99  double& tau)
100  {
101  double h,nrm_v,oneByAbsdt;
102  h = elementDiameter;
103  nrm_v=0.0;
104  for(int I=0;I<nSpace;I++)
105  nrm_v+=dH[I]*dH[I];
106  nrm_v = sqrt(nrm_v);
107  cfl = nrm_v/h;
108  oneByAbsdt = dmt;
109  //'1'
110  //tau = 1.0/(2.0*nrm_v/h + oneByAbsdt + 1.0e-8);
111  //'2'
112  tau = 1.0/sqrt(4.0*nrm_v*nrm_v/(h*h) + oneByAbsdt*oneByAbsdt + 1.0e-8);
113  }
114 
115  inline
116  void calculateSubgridError_tau( const double G[nSpace*nSpace],
117  const double Ai[nSpace],
118  double& tau_v,
119  double& q_cfl)
120  {
121  double v_d_Gv=0.0;
122  for(int I=0;I<nSpace;I++)
123  for (int J=0;J<nSpace;J++)
124  v_d_Gv += Ai[I]*G[I*nSpace+J]*Ai[J];
125 
126  tau_v = 1.0/(sqrt(v_d_Gv) + 1.0e-8);
127  }
128 
129 #undef CKDEBUG
130  void calculateResidual(arguments_dict& args,
131  bool useExact)
132  {
133  xt::pyarray<double>& mesh_trial_ref = args.array<double>("mesh_trial_ref");
134  xt::pyarray<double>& mesh_grad_trial_ref = args.array<double>("mesh_grad_trial_ref");
135  xt::pyarray<double>& mesh_dof = args.array<double>("mesh_dof");
136  xt::pyarray<int>& mesh_l2g = args.array<int>("mesh_l2g");
137  xt::pyarray<double>& x_ref = args.array<double>("x_ref");
138  xt::pyarray<double>& dV_ref = args.array<double>("dV_ref");
139  xt::pyarray<double>& u_trial_ref = args.array<double>("u_trial_ref");
140  xt::pyarray<double>& u_grad_trial_ref = args.array<double>("u_grad_trial_ref");
141  xt::pyarray<double>& u_test_ref = args.array<double>("u_test_ref");
142  xt::pyarray<double>& u_grad_test_ref = args.array<double>("u_grad_test_ref");
143  xt::pyarray<double>& mesh_trial_trace_ref = args.array<double>("mesh_trial_trace_ref");
144  xt::pyarray<double>& mesh_grad_trial_trace_ref = args.array<double>("mesh_grad_trial_trace_ref");
145  xt::pyarray<double>& dS_ref = args.array<double>("dS_ref");
146  xt::pyarray<double>& u_trial_trace_ref = args.array<double>("u_trial_trace_ref");
147  xt::pyarray<double>& u_grad_trial_trace_ref = args.array<double>("u_grad_trial_trace_ref");
148  xt::pyarray<double>& u_test_trace_ref = args.array<double>("u_test_trace_ref");
149  xt::pyarray<double>& u_grad_test_trace_ref = args.array<double>("u_grad_test_trace_ref");
150  xt::pyarray<double>& normal_ref = args.array<double>("normal_ref");
151  xt::pyarray<double>& boundaryJac_ref = args.array<double>("boundaryJac_ref");
152  int nElements_global = args.scalar<int>("nElements_global");
153  double useMetrics = args.scalar<double>("useMetrics");
154  double alphaBDF = args.scalar<double>("alphaBDF");
155  double epsFact_redist = args.scalar<double>("epsFact_redist");
156  double backgroundDiffusionFactor = args.scalar<double>("backgroundDiffusionFactor");
157  double weakDirichletFactor = args.scalar<double>("weakDirichletFactor");
158  int freezeLevelSet = args.scalar<int>("freezeLevelSet");
159  int useTimeIntegration = args.scalar<int>("useTimeIntegration");
160  int lag_shockCapturing = args.scalar<int>("lag_shockCapturing");
161  int lag_subgridError = args.scalar<int>("lag_subgridError");
162  double shockCapturingDiffusion = args.scalar<double>("shockCapturingDiffusion");
163  xt::pyarray<int>& u_l2g = args.array<int>("u_l2g");
164  xt::pyarray<double>& elementDiameter = args.array<double>("elementDiameter");
165  xt::pyarray<double>& nodeDiametersArray = args.array<double>("nodeDiametersArray");
166  xt::pyarray<double>& u_dof = args.array<double>("u_dof");
167  xt::pyarray<double>& phi_dof = args.array<double>("phi_dof");
168  xt::pyarray<double>& phi_ls = args.array<double>("phi_ls");
169  xt::pyarray<double>& q_m = args.array<double>("q_m");
170  xt::pyarray<double>& q_u = args.array<double>("q_u");
171  xt::pyarray<double>& q_n = args.array<double>("q_n");
172  xt::pyarray<double>& q_dH = args.array<double>("q_dH");
173  xt::pyarray<double>& u_weak_internal_bc_dofs = args.array<double>("u_weak_internal_bc_dofs");
174  xt::pyarray<double>& q_m_betaBDF = args.array<double>("q_m_betaBDF");
175  xt::pyarray<double>& q_dH_last = args.array<double>("q_dH_last");
176  xt::pyarray<double>& q_cfl = args.array<double>("q_cfl");
177  xt::pyarray<double>& q_numDiff_u = args.array<double>("q_numDiff_u");
178  xt::pyarray<double>& q_numDiff_u_last = args.array<double>("q_numDiff_u_last");
179  xt::pyarray<int>& weakDirichletConditionFlags = args.array<int>("weakDirichletConditionFlags");
180  int offset_u = args.scalar<int>("offset_u");
181  int stride_u = args.scalar<int>("stride_u");
182  xt::pyarray<double>& globalResidual = args.array<double>("globalResidual");
183  int nExteriorElementBoundaries_global = args.scalar<int>("nExteriorElementBoundaries_global");
184  xt::pyarray<int>& exteriorElementBoundariesArray = args.array<int>("exteriorElementBoundariesArray");
185  xt::pyarray<int>& elementBoundaryElementsArray = args.array<int>("elementBoundaryElementsArray");
186  xt::pyarray<int>& elementBoundaryLocalElementBoundariesArray = args.array<int>("elementBoundaryLocalElementBoundariesArray");
187  xt::pyarray<double>& ebqe_phi_ls_ext = args.array<double>("ebqe_phi_ls_ext");
188  xt::pyarray<int>& isDOFBoundary_u = args.array<int>("isDOFBoundary_u");
189  xt::pyarray<double>& ebqe_bc_u_ext = args.array<double>("ebqe_bc_u_ext");
190  xt::pyarray<double>& ebqe_u = args.array<double>("ebqe_u");
191  xt::pyarray<double>& ebqe_n = args.array<double>("ebqe_n");
192  int ELLIPTIC_REDISTANCING = args.scalar<int>("ELLIPTIC_REDISTANCING");
193  double backgroundDissipationEllipticRedist = args.scalar<double>("backgroundDissipationEllipticRedist");
194  xt::pyarray<double>& lumped_qx = args.array<double>("lumped_qx");
195  xt::pyarray<double>& lumped_qy = args.array<double>("lumped_qy");
196  xt::pyarray<double>& lumped_qz = args.array<double>("lumped_qz");
197  double alpha = args.scalar<double>("alpha");
198  double circbc=0.0,circ=0.0;
199  gf.useExact=useExact;
200  gfu.useExact=useExact;
201 #ifdef CKDEBUG
202  std::cout<<"stuff"<<"\t"
203  <<alphaBDF<<"\t"
204  <<epsFact_redist<<"\t"
205  <<freezeLevelSet<<"\t"
206  <<useTimeIntegration<<"\t"
207  <<lag_shockCapturing<< "\t"
208  <<lag_subgridError<<std::endl;
209 #endif
210  //
211  //loop over elements to compute volume integrals and load them into element and global residual
212  //
213  //eN is the element index
214  //eN_k is the quadrature point index for a scalar
215  //eN_k_nSpace is the quadrature point index for a vector
216  //eN_i is the element test function index
217  //eN_j is the element trial function index
218  //eN_k_j is the quadrature point index for a trial function
219  //eN_k_i is the quadrature point index for a trial function
220  double timeIntegrationScale = 1.0;
221  if (useTimeIntegration == 0)
222  timeIntegrationScale = 0.0;
223  double lag_shockCapturingScale = 1.0;
224  if (lag_shockCapturing == 0)
225  lag_shockCapturingScale = 0.0;
226  for(int eN=0;eN<nElements_global;eN++)
227  {
228  //declare local storage for element residual and initialize
229  register int dummy_l2g[nDOF_mesh_trial_element];
230  register double elementResidual_u[nDOF_test_element],element_phi[nDOF_trial_element];
231  double epsilon_redist,h_phi, dir[nSpace], norm;
232  for (int i=0;i<nDOF_test_element;i++)
233  {
234  register int eN_i=eN*nDOF_trial_element+i;
235  elementResidual_u[i]=0.0;
236  element_phi[i] = phi_dof.data()[u_l2g.data()[eN_i]];
237  dummy_l2g[i] = i;
238  }//i
239  double element_nodes[nDOF_mesh_trial_element*3];
240  for (int i=0;i<nDOF_mesh_trial_element;i++)
241  {
242  register int eN_i=eN*nDOF_mesh_trial_element+i;
243  for(int I=0;I<3;I++)
244  element_nodes[i*3 + I] = mesh_dof.data()[mesh_l2g.data()[eN_i]*3 + I];
245  }//i
246  gf.calculate(element_phi, element_nodes, x_ref.data(),false);
247  /* for (int i=0;i<nDOF_test_element;i++) */
248  /* { */
249  /* register int eN_i=eN*nDOF_trial_element+i; */
250  /* double eps=1.0e-4; */
251  /* if((fabs(gf.exact.phi_dof_corrected[i]) < eps) || */
252  /* (fabs(phi_dof.data()[u_l2g.data()[eN_i]]) < eps)) */
253  /* std::cout<<"Warning "<<gf.exact.phi_dof_corrected[i]<<'\t' */
254  /* <<phi_dof.data()[u_l2g.data()[eN_i]]<<'\t' */
255  /* <<u_l2g.data()[eN_i]<<std::endl; */
256  /* } */
257  //loop over quadrature points and compute integrands
258  for (int k=0;k<nQuadraturePoints_element;k++)
259  {
260  gf.set_quad(k);
261  //compute indeces and declare local storage
262  register int eN_k = eN*nQuadraturePoints_element+k,
263  eN_k_nSpace = eN_k*nSpace,
264  eN_nDOF_trial_element = eN*nDOF_trial_element;
265  register double u=0.0,grad_u[nSpace],u0=0.0, grad_phi[nSpace],
266  m=0.0,dm=0.0,
267  H=0.0,dH[nSpace],
268  m_t=0.0,dm_t=0.0,
269  r=0.0,
270  dH_tau[nSpace],//dH if not lagging or q_dH_last if lagging tau
271  dH_strong[nSpace],//dH if not lagging or q_dH_last if lagging strong residual and adjoint
272  pdeResidual_u=0.0,
273  Lstar_u[nDOF_test_element],
274  subgridError_u=0.0,
275  tau=0.0,tau0=0.0,tau1=0.0,
276  numDiff0=0.0,numDiff1=0.0,
277  nu_sc=0.0,
278  jac[nSpace*nSpace],
279  jacDet,
280  jacInv[nSpace*nSpace],
281  u_grad_trial[nDOF_trial_element*nSpace],
282  u_test_dV[nDOF_trial_element],
283  u_grad_test_dV[nDOF_test_element*nSpace],
284  dV,x,y,z,
285  G[nSpace*nSpace],G_dd_G,tr_G;
286  ck.calculateMapping_element(eN,
287  k,
288  mesh_dof.data(),
289  mesh_l2g.data(),
290  mesh_trial_ref.data(),
291  mesh_grad_trial_ref.data(),
292  jac,
293  jacDet,
294  jacInv,
295  x,y,z);
296  ck.calculateH_element(eN,
297  k,
298  nodeDiametersArray.data(),
299  mesh_l2g.data(),
300  mesh_trial_ref.data(),
301  h_phi);
302  //get the physical integration weight
303  dV = fabs(jacDet)*dV_ref.data()[k];
304  ck.calculateG(jacInv,G,G_dd_G,tr_G);
305 
306  //get the trial function gradients
307  ck.gradTrialFromRef(&u_grad_trial_ref.data()[k*nDOF_trial_element*nSpace],jacInv,u_grad_trial);
308  //get the solution
309  ck.valFromDOF(u_dof.data(),&u_l2g.data()[eN_nDOF_trial_element],&u_trial_ref.data()[k*nDOF_trial_element],u);
310  ck.valFromDOF(gf.exact.phi_dof_corrected,dummy_l2g,&u_trial_ref.data()[k*nDOF_trial_element],u0);
311  if (freezeLevelSet)
312  {
313  u0 = phi_ls.data()[eN_k];
314  }
315  /* double DX=(x-0.5),DY=(y-0.75); */
316  /* double radius = std::sqrt(DX*DX+DY*DY); */
317  /* double theta = std::atan2(DY,DX); */
318  /* double kp=10.0, scaling=1.0, rp=1; */
319  /* u0 = scaling*std::pow((0.15+(0.015/2.)*std::cos(kp*theta) - radius),rp); */
320 
321  //u0 = 0.15 - std::sqrt(DX*DX + DY*DY);
322  //u0 = phi_ls.data()[eN_k];//cek debug--set to exact input
323  /* //cek debug */
324  /* double u0_test; */
325  /* ck.valFromDOF(phi_dof.data(),&u_l2g.data()[eN_nDOF_trial_element],&u_trial_ref.data()[k*nDOF_trial_element],u0_test); */
326  /* if (u0 != u0_test) */
327  /* std::cout<<"eN "<<eN<<" k "<<k<<" u0 "<<u0<<" u0_test "<<u0_test<<std::endl; */
328  //get the solution gradients
329  ck.gradFromDOF(u_dof.data(),&u_l2g.data()[eN_nDOF_trial_element],u_grad_trial,grad_u);
330  //precalculate test function products with integration weights
331  for (int j=0;j<nDOF_trial_element;j++)
332  {
333  u_test_dV[j] = u_test_ref.data()[k*nDOF_trial_element+j]*dV;
334  for (int I=0;I<nSpace;I++)
335  {
336  u_grad_test_dV[j*nSpace+I] = u_grad_trial[j*nSpace+I]*dV;//cek warning won't work for Petrov-Galerkin
337  }
338  }
339  //
340  //calculate pde coefficients at quadrature points
341  //
342  /* norm = 1.0e-8; */
343  /* for (int I=0;I<nSpace;I++) */
344  /* norm += grad_u[I]*grad_u[I]; */
345  /* norm = sqrt(norm); */
346 
347  /* for (int I=0;I<nSpace;I++) */
348  /* dir[I] = grad_u[I]/norm; */
349 
350  /* ck.calculateGScale(G,dir,h_phi); */
351 
352  epsilon_redist = epsFact_redist*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter.data()[eN]);
353 
354  evaluateCoefficients(epsilon_redist,
355  u0,
356  u,
357  grad_u,
358  m,
359  dm,
360  H,
361  dH,
362  r);
363  //TODO allow not lagging of subgrid error etc,
364  //remove conditional?
365  //default no lagging
366  for (int I=0; I < nSpace; I++)
367  {
368  dH_tau[I] = dH[I];
369  dH_strong[I] = dH[I];
370  }
371  if (lag_subgridError > 0)
372  {
373  for (int I=0; I < nSpace; I++)
374  {
375  dH_tau[I] = q_dH_last.data()[eN_k_nSpace+I];
376  }
377  }
378  if (lag_subgridError > 1)
379  {
380  for (int I=0; I < nSpace; I++)
381  {
382  dH_strong[I] = q_dH_last.data()[eN_k_nSpace+I];
383  }
384  }
385  //save mass for time history and dH for subgrid error
386  //save solution for other models
387  //
388  q_m.data()[eN_k] = m;
389  q_u.data()[eN_k] = u;
390  for (int I=0;I<nSpace;I++)
391  q_n.data()[eN_k_nSpace+I] = dir[I];
392 
393  for (int I=0;I<nSpace;I++)
394  {
395  int eN_k_nSpace_I = eN_k_nSpace+I;
396  q_dH.data()[eN_k_nSpace_I] = dH[I];
397  }
398 
399  //
400  //moving mesh
401  //
402  //omit for now
403  //
404  //calculate time derivative at quadrature points
405  //
406  ck.bdf(alphaBDF,
407  q_m_betaBDF.data()[eN_k],
408  m,
409  dm,
410  m_t,
411  dm_t);
412 
413  //TODO add option to skip if not doing time integration (Newton stead-state solve)
414  m *= timeIntegrationScale; dm *= timeIntegrationScale; m_t *= timeIntegrationScale;
415  dm_t *= timeIntegrationScale;
416 #ifdef CKDEBUG
417  std::cout<<"alpha "<<alphaBDF<<"\t"<<q_m_betaBDF.data()[eN_k]<<"\t"<<m_t<<'\t'<<m<<'\t'<<alphaBDF*m<<std::endl;
418 #endif
419  //
420  //calculate subgrid error (strong residual and adjoint)
421  //
422  //calculate strong residual
423  pdeResidual_u = ck.Mass_strong(m_t) +
424  ck.Hamiltonian_strong(dH_strong,grad_u) + //would need dH if not lagging
425  ck.Reaction_strong(r);
426 #ifdef CKDEBUG
427  std::cout<<"dH_strong "<<dH_strong[0]<<'\t'<<dH_strong[1]<<'\t'<<dH_strong[2]<<std::endl;
428 #endif
429  //calculate adjoint
430  for (int i=0;i<nDOF_test_element;i++)
431  {
432  //register int eN_k_i_nSpace = (eN_k*nDOF_trial_element+i)*nSpace;
433  register int i_nSpace=i*nSpace;
434  Lstar_u[i] = ck.Hamiltonian_adjoint(dH_strong,&u_grad_test_dV[i_nSpace]);
435  //reaction is constant
436  }
437  //calculate tau and tau*Res
438  calculateSubgridError_tau(elementDiameter.data()[eN],
439  dm_t,dH_tau,
440  q_cfl.data()[eN_k],
441  tau0);
442  calculateSubgridError_tau(G,
443  dH_tau,
444  tau1,
445  q_cfl.data()[eN_k]);
446 
447  tau = useMetrics*tau1+(1.0-useMetrics)*tau0;
448 
449  //std::cout<<tau<<std::endl;
450 
451 
452  subgridError_u = -tau*pdeResidual_u;
453  //
454  //calculate shock capturing diffusion
455  //
456  ck.calculateNumericalDiffusion(shockCapturingDiffusion,elementDiameter.data()[eN],pdeResidual_u,grad_u,numDiff0);
457  ck.calculateNumericalDiffusion(shockCapturingDiffusion,G,pdeResidual_u,grad_u,numDiff1);
458 
459  q_numDiff_u.data()[eN_k] = useMetrics*numDiff1+(1.0-useMetrics)*numDiff0;
460 
461 #ifdef CKDEBUG
462  std::cout<<"q_numDiff_u[eN_k] "<<q_numDiff_u.data()[eN_k]<<" q_numDiff_u_last[eN_k] "<<q_numDiff_u_last.data()[eN_k]<<" lag "<<lag_shockCapturingScale<<std::endl;
463 #endif
464  nu_sc = q_numDiff_u.data()[eN_k]*(1.0-lag_shockCapturingScale) + q_numDiff_u_last.data()[eN_k]*lag_shockCapturingScale;
465  /* double epsilon_background_diffusion = 3.0*h_phi;//2.0*epsFact_redist*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter.data()[eN]); */
466  /* if (fabs(u0) > epsilon_background_diffusion) */
467  /* nu_sc += backgroundDiffusionFactor*h_phi; */
468 
469  double epsilon_background_diffusion = 2.0*epsFact_redist*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter.data()[eN]);
470  if (fabs(phi_ls.data()[eN_k]) > epsilon_background_diffusion)
471  nu_sc += backgroundDiffusionFactor*elementDiameter.data()[eN]; //
472  //update element residual
473  //
474  for(int i=0;i<nDOF_test_element;i++)
475  {
476  register int i_nSpace = i*nSpace;
477  double FREEZE=double(freezeLevelSet);
478  //assert(FREEZE==0.0);
479  //register int eN_k_i=eN_k*nDOF_test_element+i;
480  //register int eN_k_i_nSpace = eN_k_i*nSpace;
481 
482 #ifdef CKDEBUG
483  std::cout<<"shock capturing input nu_sc "<<nu_sc<<'\t'<<grad_u[0]<<'\t'<<grad_u[1]<<'\t'<<grad_u[1]<<'\t'<<u_grad_test_dV[i_nSpace]<<std::endl;
484 #endif
485  //std::cout<<element_phi[i]<<'\t'<<element_nodes[i*3+0]<<'\t'<<element_nodes[i*3+1]<<std::endl;
486  //std::cout<<"eN_k "<<eN_k<<" D "<<gf.D(epsilon_redist,phi_ls.data()[eN_k])<<std::endl;
487  circbc += ck.Reaction_weak(gf.D(epsilon_redist,u0)*(u-u0), u_test_dV[i]);
488  elementResidual_u[i] += ck.Mass_weak(m_t,u_test_dV[i]) +
489  ck.Hamiltonian_weak(H,u_test_dV[i]) +
490  ck.Reaction_weak(r,u_test_dV[i]) +
491  (1.0-FREEZE)*(weakDirichletFactor/h_phi)*ck.Reaction_weak(gf.D(epsilon_redist,u0)*(u0-u),
492  u_test_dV[i]) +
493  ck.SubgridError(subgridError_u,Lstar_u[i]) +
494  ck.NumericalDiffusion(nu_sc,grad_u,&u_grad_test_dV[i_nSpace]);
495 #ifdef CKDEBUG
496  std::cout<<ck.Mass_weak(m_t,u_test_dV[i])<<'\t'
497  <<ck.Hamiltonian_weak(H,u_test_dV[i]) <<'\t'
498  <<ck.Reaction_weak(r,u_test_dV[i])<<'\t'
499  <<ck.SubgridError(subgridError_u,Lstar_u[i])<<'\t'
500  <<ck.NumericalDiffusion(nu_sc,grad_u,&u_grad_test_dV[i_nSpace])<<std::endl;
501 #endif
502  }//i
503  //
504  }//k
505  //
506  //apply weak constraints for unknowns near zero level set
507  //
508  //
509  if (freezeLevelSet)
510  {
511  for (int j = 0; j < nDOF_trial_element; j++)
512  {
513  const int eN_j = eN*nDOF_trial_element+j;
514  const int J = u_l2g.data()[eN_j];
515  //if (weakDirichletConditionFlags.data()[J] == 1)
516  if (fabs(u_weak_internal_bc_dofs.data()[J]) < epsilon_redist)
517  {
518  elementResidual_u[j] = (u_dof.data()[J]-u_weak_internal_bc_dofs.data()[J])*weakDirichletFactor*elementDiameter.data()[eN];
519  }
520  }//j
521  }//freeze
522 
523  //
524  //load element into global residual and save element residual
525  //
526  for(int i=0;i<nDOF_test_element;i++)
527  {
528  register int eN_i=eN*nDOF_test_element+i;
529 
530 #ifdef CKDEBUG
531  std::cout<<"element residual i = "<<i<<"\t"<<elementResidual_u[i]<<std::endl;
532 #endif
533  globalResidual.data()[offset_u+stride_u*u_l2g.data()[eN_i]]+=elementResidual_u[i];
534  }//i
535  }//elements
536  //
537  //loop over exterior element boundaries
538  //
539  //ebNE is the Exterior element boundary INdex
540  //ebN is the element boundary INdex
541  //eN is the element index
542  for (int ebNE = 0; ebNE < nExteriorElementBoundaries_global; ebNE++)
543  {
544  register int ebN = exteriorElementBoundariesArray.data()[ebNE],
545  eN = elementBoundaryElementsArray.data()[ebN*2+0],
546  ebN_local = elementBoundaryLocalElementBoundariesArray.data()[ebN*2+0],
547  eN_nDOF_trial_element = eN*nDOF_trial_element;
548  double epsilon_redist, h_phi;
549  for (int kb=0;kb<nQuadraturePoints_elementBoundary;kb++)
550  {
551  register int ebNE_kb = ebNE*nQuadraturePoints_elementBoundary+kb,
552  ebNE_kb_nSpace = ebNE_kb*nSpace,
553  ebN_local_kb = ebN_local*nQuadraturePoints_elementBoundary+kb,
554  ebN_local_kb_nSpace = ebN_local_kb*nSpace;
555  register double u_ext=0.0,
556  grad_u_ext[nSpace],
557  jac_ext[nSpace*nSpace],
558  jacDet_ext,
559  jacInv_ext[nSpace*nSpace],
560  boundaryJac[nSpace*(nSpace-1)],
561  metricTensor[(nSpace-1)*(nSpace-1)],
562  metricTensorDetSqrt,
563  u_test_dS[nDOF_test_element],
564  u_grad_trial_trace[nDOF_trial_element*nSpace],
565  normal[nSpace],x_ext,y_ext,z_ext,
566  dir[nSpace],norm;
567  ck.calculateMapping_elementBoundary(eN,
568  ebN_local,
569  kb,
570  ebN_local_kb,
571  mesh_dof.data(),
572  mesh_l2g.data(),
573  mesh_trial_trace_ref.data(),
574  mesh_grad_trial_trace_ref.data(),
575  boundaryJac_ref.data(),
576  jac_ext,
577  jacDet_ext,
578  jacInv_ext,
579  boundaryJac,
580  metricTensor,
581  metricTensorDetSqrt,
582  normal_ref.data(),
583  normal,
584  x_ext,y_ext,z_ext);
585  //compute shape and solution information
586  //shape
587  ck.gradTrialFromRef(&u_grad_trial_trace_ref.data()[ebN_local_kb_nSpace*nDOF_trial_element],jacInv_ext,u_grad_trial_trace);
588  //solution and gradients
589  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);
590  ck.gradFromDOF(u_dof.data(),&u_l2g.data()[eN_nDOF_trial_element],u_grad_trial_trace,grad_u_ext);
591  norm = 1.0e-8;
592  for (int I=0;I<nSpace;I++)
593  norm += grad_u_ext[I]*grad_u_ext[I];
594  norm = sqrt(norm);
595  for (int I=0;I<nSpace;I++)
596  dir[I] = grad_u_ext[I]/norm;
597  //save for other models
598  ebqe_u.data()[ebNE_kb] = u_ext;
599  for (int I=0;I<nSpace;I++)
600  ebqe_n.data()[ebNE_kb_nSpace+I] = dir[I];
601  }//kb
602  }//ebNE
603  //std::cout<<"Circ Res"<<circbc<<'\t'<<circ<<std::endl;
604  }
605 
606  //for now assumes that using time integration
607  //and so lags stabilization and subgrid error
608  //extern "C" void calculateJacobian_RDLSV2(int nElements_global,
609  void calculateJacobian(arguments_dict& args,
610  bool useExact)
611  {
612  xt::pyarray<double>& mesh_trial_ref = args.array<double>("mesh_trial_ref");
613  xt::pyarray<double>& mesh_grad_trial_ref = args.array<double>("mesh_grad_trial_ref");
614  xt::pyarray<double>& mesh_dof = args.array<double>("mesh_dof");
615  xt::pyarray<int>& mesh_l2g = args.array<int>("mesh_l2g");
616  xt::pyarray<double>& x_ref = args.array<double>("x_ref");
617  xt::pyarray<double>& dV_ref = args.array<double>("dV_ref");
618  xt::pyarray<double>& u_trial_ref = args.array<double>("u_trial_ref");
619  xt::pyarray<double>& u_grad_trial_ref = args.array<double>("u_grad_trial_ref");
620  xt::pyarray<double>& u_test_ref = args.array<double>("u_test_ref");
621  xt::pyarray<double>& u_grad_test_ref = args.array<double>("u_grad_test_ref");
622  xt::pyarray<double>& mesh_trial_trace_ref = args.array<double>("mesh_trial_trace_ref");
623  xt::pyarray<double>& mesh_grad_trial_trace_ref = args.array<double>("mesh_grad_trial_trace_ref");
624  xt::pyarray<double>& dS_ref = args.array<double>("dS_ref");
625  xt::pyarray<double>& u_trial_trace_ref = args.array<double>("u_trial_trace_ref");
626  xt::pyarray<double>& u_grad_trial_trace_ref = args.array<double>("u_grad_trial_trace_ref");
627  xt::pyarray<double>& u_test_trace_ref = args.array<double>("u_test_trace_ref");
628  xt::pyarray<double>& u_grad_test_trace_ref = args.array<double>("u_grad_test_trace_ref");
629  xt::pyarray<double>& normal_ref = args.array<double>("normal_ref");
630  xt::pyarray<double>& boundaryJac_ref = args.array<double>("boundaryJac_ref");
631  int nElements_global = args.scalar<int>("nElements_global");
632  double useMetrics = args.scalar<double>("useMetrics");
633  double alphaBDF = args.scalar<double>("alphaBDF");
634  double epsFact_redist = args.scalar<double>("epsFact_redist");
635  double backgroundDiffusionFactor = args.scalar<double>("backgroundDiffusionFactor");
636  double weakDirichletFactor = args.scalar<double>("weakDirichletFactor");
637  int freezeLevelSet = args.scalar<int>("freezeLevelSet");
638  int useTimeIntegration = args.scalar<int>("useTimeIntegration");
639  int lag_shockCapturing = args.scalar<int>("lag_shockCapturing");
640  int lag_subgridError = args.scalar<int>("lag_subgridError");
641  double shockCapturingDiffusion = args.scalar<double>("shockCapturingDiffusion");
642  xt::pyarray<int>& u_l2g = args.array<int>("u_l2g");
643  xt::pyarray<double>& elementDiameter = args.array<double>("elementDiameter");
644  xt::pyarray<double>& nodeDiametersArray = args.array<double>("nodeDiametersArray");
645  xt::pyarray<double>& u_dof = args.array<double>("u_dof");
646  xt::pyarray<double>& phi_dof = args.array<double>("phi_dof");
647  xt::pyarray<double>& u_weak_internal_bc_dofs = args.array<double>("u_weak_internal_bc_dofs");
648  xt::pyarray<double>& phi_ls = args.array<double>("phi_ls");
649  xt::pyarray<double>& q_m_betaBDF = args.array<double>("q_m_betaBDF");
650  xt::pyarray<double>& q_dH_last = args.array<double>("q_dH_last");
651  xt::pyarray<double>& q_cfl = args.array<double>("q_cfl");
652  xt::pyarray<double>& q_numDiff_u = args.array<double>("q_numDiff_u");
653  xt::pyarray<double>& q_numDiff_u_last = args.array<double>("q_numDiff_u_last");
654  xt::pyarray<int>& weakDirichletConditionFlags = args.array<int>("weakDirichletConditionFlags");
655  xt::pyarray<int>& csrRowIndeces_u_u = args.array<int>("csrRowIndeces_u_u");
656  xt::pyarray<int>& csrColumnOffsets_u_u = args.array<int>("csrColumnOffsets_u_u");
657  xt::pyarray<double>& globalJacobian = args.array<double>("globalJacobian");
658  int nExteriorElementBoundaries_global = args.scalar<int>("nExteriorElementBoundaries_global");
659  xt::pyarray<int>& exteriorElementBoundariesArray = args.array<int>("exteriorElementBoundariesArray");
660  xt::pyarray<int>& elementBoundaryElementsArray = args.array<int>("elementBoundaryElementsArray");
661  xt::pyarray<int>& elementBoundaryLocalElementBoundariesArray = args.array<int>("elementBoundaryLocalElementBoundariesArray");
662  xt::pyarray<double>& ebqe_phi_ls_ext = args.array<double>("ebqe_phi_ls_ext");
663  xt::pyarray<int>& isDOFBoundary_u = args.array<int>("isDOFBoundary_u");
664  xt::pyarray<double>& ebqe_bc_u_ext = args.array<double>("ebqe_bc_u_ext");
665  xt::pyarray<int>& csrColumnOffsets_eb_u_u = args.array<int>("csrColumnOffsets_eb_u_u");
666  int ELLIPTIC_REDISTANCING = args.scalar<int>("ELLIPTIC_REDISTANCING");
667  double backgroundDissipationEllipticRedist = args.scalar<double>("backgroundDissipationEllipticRedist");
668  double alpha = args.scalar<double>("alpha");
669  double circ=0.0;
670  gf.useExact=useExact;
671  //
672  //loop over elements to compute volume integrals and load them into the element Jacobians and global Jacobian
673  //
674  double timeIntegrationScale = 1.0;
675  if (useTimeIntegration == 0)
676  timeIntegrationScale = 0.0;
677  double lag_shockCapturingScale = 1.0;
678  if (lag_shockCapturing == 0)
679  lag_shockCapturingScale = 0.0;
680  for(int eN=0;eN<nElements_global;eN++)
681  {
682  register int dummy_l2g[nDOF_mesh_trial_element];
683  register double elementJacobian_u_u[nDOF_test_element][nDOF_trial_element],element_phi[nDOF_trial_element];
684  double epsilon_redist,h_phi, dir[nSpace], norm;
685  for (int i=0;i<nDOF_test_element;i++)
686  {
687  register int eN_i=eN*nDOF_trial_element+i;
688  element_phi[i] = phi_dof.data()[u_l2g.data()[eN_i]];
689  dummy_l2g[i] = i;
690  for (int j=0;j<nDOF_trial_element;j++)
691  {
692  elementJacobian_u_u[i][j]=0.0;
693  }
694  }
695  double element_nodes[nDOF_mesh_trial_element*3];
696  for (int i=0;i<nDOF_mesh_trial_element;i++)
697  {
698  register int eN_i=eN*nDOF_mesh_trial_element+i;
699  for(int I=0;I<3;I++)
700  element_nodes[i*3 + I] = mesh_dof.data()[mesh_l2g.data()[eN_i]*3 + I];
701  }//i
702  gf.calculate(element_phi, element_nodes, x_ref.data(),false);
703  for (int k=0;k<nQuadraturePoints_element;k++)
704  {
705  gf.set_quad(k);
706  int eN_k = eN*nQuadraturePoints_element+k, //index to a scalar at a quadrature point
707  eN_k_nSpace = eN_k*nSpace,
708  eN_nDOF_trial_element = eN*nDOF_trial_element; //index to a vector at a quadrature point
709 
710  //declare local storage
711  register double u=0.0,u0=0.0,
712  grad_u[nSpace],
713  m=0.0,dm=0.0,
714  H=0.0,dH[nSpace],
715  m_t=0.0,dm_t=0.0,r=0.0,
716  dH_tau[nSpace],//dH or dH_last if lagging for tau formula
717  dH_strong[nSpace],//dH or dH_last if lagging for strong residual and adjoint
718  dpdeResidual_u_u[nDOF_trial_element],
719  Lstar_u[nDOF_test_element],
720  dsubgridError_u_u[nDOF_trial_element],
721  tau=0.0,tau0=0.0,tau1=0.0,
722  nu_sc=0.0,
723  jac[nSpace*nSpace],
724  jacDet,
725  jacInv[nSpace*nSpace],
726  u_grad_trial[nDOF_trial_element*nSpace],
727  dV,
728  u_test_dV[nDOF_test_element],
729  u_grad_test_dV[nDOF_test_element*nSpace],
730  x,y,z,
731  G[nSpace*nSpace],G_dd_G,tr_G;
732  //
733  //calculate solution and gradients at quadrature points
734  //
735  ck.calculateMapping_element(eN,
736  k,
737  mesh_dof.data(),
738  mesh_l2g.data(),
739  mesh_trial_ref.data(),
740  mesh_grad_trial_ref.data(),
741  jac,
742  jacDet,
743  jacInv,
744  x,y,z);
745  ck.calculateH_element(eN,
746  k,
747  nodeDiametersArray.data(),
748  mesh_l2g.data(),
749  mesh_trial_ref.data(),
750  h_phi);
751  //get the physical integration weight
752  dV = fabs(jacDet)*dV_ref.data()[k];
753  ck.calculateG(jacInv,G,G_dd_G,tr_G);
754  //get the trial function gradients
755  ck.gradTrialFromRef(&u_grad_trial_ref.data()[k*nDOF_trial_element*nSpace],jacInv,u_grad_trial);
756  //get the solution
757  ck.valFromDOF(u_dof.data(),&u_l2g.data()[eN_nDOF_trial_element],&u_trial_ref.data()[k*nDOF_trial_element],u);
758  ck.valFromDOF(gf.exact.phi_dof_corrected,dummy_l2g,&u_trial_ref.data()[k*nDOF_trial_element],u0);
759  if (freezeLevelSet)
760  u0 = phi_ls.data()[eN_k];
761  //u0 = phi_ls.data()[eN_k];//cek debug
762  /* double DX=(x-0.5),DY=(y-0.75); */
763  /* double radius = std::sqrt(DX*DX+DY*DY); */
764  /* double theta = std::atan2(DY,DX); */
765  /* double kp=10.0, scaling=1.0, rp=1; */
766  /* u0 = scaling*std::pow((0.15+(0.015/2.)*std::cos(kp*theta) - radius),rp); */
767  //u0 = 0.15 - std::sqrt(DX*DX + DY*DY);
768  //u0 = phi_ls.data()[eN_k];//cek debug--set to exact input
769  /* double u0_test=0; */
770  /* ck.valFromDOF(phi_dof.data(),&u_l2g.data()[eN_nDOF_trial_element],&u_trial_ref.data()[k*nDOF_trial_element],u0_test); */
771  /* if (u0 != u0_test) */
772  /* std::cout<<"JAC eN "<<eN<<" k "<<k<<" u0 "<<u0<<" u0_test "<<u0_test<<std::endl; */
773 
774  //get the solution gradients
775  ck.gradFromDOF(u_dof.data(),&u_l2g.data()[eN_nDOF_trial_element],u_grad_trial,grad_u);
776  //precalculate test function products with integration weights
777  for (int j=0;j<nDOF_trial_element;j++)
778  {
779  u_test_dV[j] = u_test_ref.data()[k*nDOF_trial_element+j]*dV;
780  for (int I=0;I<nSpace;I++)
781  {
782  u_grad_test_dV[j*nSpace+I] = u_grad_trial[j*nSpace+I]*dV;//cek warning won't work for Petrov-Galerkin
783  }
784  }
785  //
786  //calculate pde coefficients and derivatives at quadrature points
787  //
788  /* norm = 1.0e-8; */
789  /* for (int I=0;I<nSpace;I++) */
790  /* norm += grad_u[I]*grad_u[I]; */
791  /* norm = sqrt(norm); */
792  /* for (int I=0;I<nSpace;I++) */
793  /* dir[I] = grad_u[I]/norm; */
794 
795  /* ck.calculateGScale(G,dir,h_phi); */
796 
797  epsilon_redist = epsFact_redist*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter.data()[eN]);
798 
799  evaluateCoefficients(epsilon_redist,
800  u0,
801  u,
802  grad_u,
803  m,
804  dm,
805  H,
806  dH,
807  r);
808  //TODO allow not lagging of subgrid error etc
809  //remove conditional?
810  //default no lagging
811  for (int I=0; I < nSpace; I++)
812  {
813  dH_tau[I] = dH[I];
814  dH_strong[I] = dH[I];
815  }
816  if (lag_subgridError > 0)
817  {
818  for (int I=0; I < nSpace; I++)
819  {
820  dH_tau[I] = q_dH_last.data()[eN_k_nSpace+I];
821  }
822  }
823  if (lag_subgridError > 1)
824  {
825  for (int I=0; I < nSpace; I++)
826  {
827  dH_strong[I] = q_dH_last.data()[eN_k_nSpace+I];
828  }
829  }
830  //
831  //moving mesh
832  //
833  //omit for now
834  //
835  //calculate time derivatives
836  //
837  ck.bdf(alphaBDF,
838  q_m_betaBDF.data()[eN_k],
839  m,
840  dm,
841  m_t,
842  dm_t);
843  //TODO add option to skip if not doing time integration (Newton stead-state solve)
844  m *= timeIntegrationScale; dm *= timeIntegrationScale; m_t *= timeIntegrationScale;
845  dm_t *= timeIntegrationScale;
846 
847  //
848  //calculate subgrid error contribution to the Jacobian (strong residual, adjoint, jacobian of strong residual)
849  //
850  //calculate the adjoint times the test functions
851  for (int i=0;i<nDOF_test_element;i++)
852  {
853  //int eN_k_i_nSpace = (eN_k*nDOF_trial_element+i)*nSpace;
854  int i_nSpace=i*nSpace;
855 
856  Lstar_u[i]=ck.Hamiltonian_adjoint(dH_strong,&u_grad_test_dV[i_nSpace]);
857 
858  }
859  //calculate the Jacobian of strong residual
860  for (int j=0;j<nDOF_trial_element;j++)
861  {
862  //int eN_k_j=eN_k*nDOF_trial_element+j;
863  //int eN_k_j_nSpace = eN_k_j*nSpace;
864  int j_nSpace = j*nSpace;
865  dpdeResidual_u_u[j]=ck.MassJacobian_strong(dm_t,u_trial_ref.data()[k*nDOF_trial_element+j]) +
866  ck.HamiltonianJacobian_strong(dH_strong,&u_grad_trial[j_nSpace]);
867  }
868  //tau and tau*Res
869  calculateSubgridError_tau(elementDiameter.data()[eN],
870  dm_t,
871  dH_tau,
872  q_cfl.data()[eN_k],
873  tau0);
874  calculateSubgridError_tau(G,
875  dH_tau,
876  tau1,
877  q_cfl.data()[eN_k]);
878 
879  tau = useMetrics*tau1+(1.0-useMetrics)*tau0;
880 
881  for (int j=0;j<nDOF_trial_element;j++)
882  dsubgridError_u_u[j] = -tau*dpdeResidual_u_u[j];
883 
884  nu_sc = q_numDiff_u.data()[eN_k]*(1.0-lag_shockCapturingScale) + q_numDiff_u_last.data()[eN_k]*lag_shockCapturingScale;
885  /* double epsilon_background_diffusion = 3.0*epsFact_redist*h_phi; */
886  /* if (fabs(u0) > epsilon_background_diffusion) */
887  /* nu_sc += backgroundDiffusionFactor*h_phi; */
888  double epsilon_background_diffusion = 2.0*epsFact_redist*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter.data()[eN]);
889  if (fabs(phi_ls.data()[eN_k]) > epsilon_background_diffusion)
890  nu_sc += backgroundDiffusionFactor*elementDiameter.data()[eN]; //
891  for(int i=0;i<nDOF_test_element;i++)
892  {
893  //int eN_k_i=eN_k*nDOF_test_element+i;
894  //int eN_k_i_nSpace=eN_k_i*nSpace;
895  circ += ck.Reaction_weak(gf.D(epsilon_redist, phi_ls.data()[eN_k]), u_test_dV[i]);
896  for(int j=0;j<nDOF_trial_element;j++)
897  {
898  //int eN_k_j=eN_k*nDOF_trial_element+j;
899  //int eN_k_j_nSpace = eN_k_j*nSpace;
900  int j_nSpace = j*nSpace;
901  int i_nSpace = i*nSpace;
902  double FREEZE=double(freezeLevelSet);
903  //std::cout<<element_phi[i]<<'\t'<<element_nodes[i*3+0]<<'\t'<<element_nodes[i*3+1]<<std::endl;
904  //std::cout<<"eN_k "<<eN_k<<" D-J "<<gf.D(epsilon_redist,phi_ls.data()[eN_k])<<std::endl;
905  elementJacobian_u_u[i][j] += ck.MassJacobian_weak(dm_t,u_trial_ref.data()[k*nDOF_trial_element+j],u_test_dV[i]) +
906  ck.HamiltonianJacobian_weak(dH,&u_grad_trial[j_nSpace],u_test_dV[i]) +
907  (1.0-FREEZE)*(weakDirichletFactor/h_phi)*ck.ReactionJacobian_weak(-gf.D(epsilon_redist,u0),
908  u_trial_ref.data()[k*nDOF_trial_element+j],
909  u_test_dV[i]) +
910  ck.SubgridErrorJacobian(dsubgridError_u_u[j],Lstar_u[i]) +
911  ck.NumericalDiffusionJacobian(nu_sc,&u_grad_trial[j_nSpace],&u_grad_test_dV[i_nSpace]);
912  }//j
913  }//i
914  }//k
915  //
916  //load into element Jacobian into global Jacobian
917  //
918 
919  //now try to account for weak dirichlet conditions in interior (frozen level set values)
920  if (freezeLevelSet)
921  {
922  //assume correspondence between dof and equations
923  for (int j = 0; j < nDOF_trial_element; j++)
924  {
925  const int J = u_l2g.data()[eN*nDOF_trial_element+j];
926  if (fabs(u_weak_internal_bc_dofs.data()[J]) < epsilon_redist)
927  //if (weakDirichletConditionFlags[J] == 1)
928  //if (fabs(gf.exact.phi_dof_corrected[j]) < epsilon_redist)
929  {
930  for (int jj=0; jj < nDOF_trial_element; jj++)
931  elementJacobian_u_u[j][jj] = 0.0;
932  elementJacobian_u_u[j][j] = weakDirichletFactor*elementDiameter.data()[eN];
933  }
934  }
935  }
936  for (int i=0;i<nDOF_test_element;i++)
937  {
938  int eN_i = eN*nDOF_test_element+i;
939  for (int j=0;j<nDOF_trial_element;j++)
940  {
941  int eN_i_j = eN_i*nDOF_trial_element+j;
942 #ifdef CKDEBUG
943  std::cout<<"element jacobian i = "<<i<<"\t"<<"j = "<<j<<"\t"<<elementJacobian_u_u[i][j]<<std::endl;
944 #endif
945  globalJacobian.data()[csrRowIndeces_u_u.data()[eN_i] + csrColumnOffsets_u_u.data()[eN_i_j]] += elementJacobian_u_u[i][j];
946  }//j
947  }//i
948  }//elements
949  /* //cek todo should get rid of this, see res */
950  /* // */
951  /* //loop over exterior element boundaries to compute the surface integrals and load them into the global Jacobian */
952  /* // */
953  /* gf.useExact=false;//exact Heaviside integration not implemented for boundaries yet */
954  /* for (int ebNE = 0; ebNE < nExteriorElementBoundaries_global; ebNE++) */
955  /* { */
956  /* register int ebN = exteriorElementBoundariesArray.data()[ebNE]; */
957  /* register int eN = elementBoundaryElementsArray.data()[ebN*2+0], */
958  /* ebN_local = elementBoundaryLocalElementBoundariesArray.data()[ebN*2+0], */
959  /* eN_nDOF_trial_element = eN*nDOF_trial_element; */
960  /* double epsilon_redist,h_phi; */
961  /* for (int kb=0;kb<nQuadraturePoints_elementBoundary;kb++) */
962  /* { */
963  /* register int ebNE_kb = ebNE*nQuadraturePoints_elementBoundary+kb; */
964  /* //register int ebNE_kb_nSpace = ebNE_kb*nSpace; */
965  /* register int ebN_local_kb = ebN_local*nQuadraturePoints_elementBoundary+kb; */
966  /* register int ebN_local_kb_nSpace = ebN_local_kb*nSpace; */
967 
968  /* register double u_ext=0.0, */
969  /* grad_u_ext[nSpace], */
970  /* m_ext=0.0, */
971  /* dm_ext=0.0, */
972  /* H_ext=0.0, */
973  /* dH_ext[nSpace], */
974  /* r_ext=0.0, */
975  /* // flux_ext=0.0, */
976  /* dflux_u_u_ext=0.0, */
977  /* bc_u_ext=0.0, */
978  /* bc_grad_u_ext[nSpace], */
979  /* bc_m_ext=0.0, */
980  /* bc_dm_ext=0.0, */
981  /* bc_H_ext=0.0, */
982  /* bc_dH_ext[nSpace], */
983  /* bc_r_ext=0.0, */
984  /* fluxJacobian_u_u[nDOF_trial_element], */
985  /* jac_ext[nSpace*nSpace], */
986  /* jacDet_ext, */
987  /* jacInv_ext[nSpace*nSpace], */
988  /* boundaryJac[nSpace*(nSpace-1)], */
989  /* metricTensor[(nSpace-1)*(nSpace-1)], */
990  /* metricTensorDetSqrt, */
991  /* dS, */
992  /* u_test_dS[nDOF_test_element], */
993  /* u_grad_trial_trace[nDOF_trial_element*nSpace], */
994  /* normal[nSpace],x_ext,y_ext,z_ext, */
995  /* G[nSpace*nSpace],G_dd_G,tr_G, dir[nSpace],norm; */
996  /* // */
997  /* //calculate the solution and gradients at quadrature points */
998  /* // */
999  /* ck.calculateMapping_elementBoundary(eN, */
1000  /* ebN_local, */
1001  /* kb, */
1002  /* ebN_local_kb, */
1003  /* mesh_dof.data(), */
1004  /* mesh_l2g.data(), */
1005  /* mesh_trial_trace_ref.data(), */
1006  /* mesh_grad_trial_trace_ref.data(), */
1007  /* boundaryJac_ref.data(), */
1008  /* jac_ext, */
1009  /* jacDet_ext, */
1010  /* jacInv_ext, */
1011  /* boundaryJac, */
1012  /* metricTensor, */
1013  /* metricTensorDetSqrt, */
1014  /* normal_ref.data(), */
1015  /* normal, */
1016  /* x_ext,y_ext,z_ext); */
1017  /* dS = metricTensorDetSqrt*dS_ref.data()[kb]; */
1018  /* ck.calculateG(jacInv_ext,G,G_dd_G,tr_G); */
1019  /* //compute shape and solution information */
1020  /* //shape */
1021  /* ck.gradTrialFromRef(&u_grad_trial_trace_ref.data()[ebN_local_kb_nSpace*nDOF_trial_element],jacInv_ext,u_grad_trial_trace); */
1022  /* //solution and gradients */
1023  /* 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); */
1024  /* ck.gradFromDOF(u_dof.data(),&u_l2g.data()[eN_nDOF_trial_element],u_grad_trial_trace,grad_u_ext); */
1025  /* //precalculate test function products with integration weights */
1026  /* for (int j=0;j<nDOF_trial_element;j++) */
1027  /* { */
1028  /* u_test_dS[j] = u_test_trace_ref.data()[ebN_local_kb*nDOF_test_element+j]*dS; */
1029  /* } */
1030 
1031  /* norm = 1.0e-8; */
1032  /* for (int I=0;I<nSpace;I++) */
1033  /* norm += grad_u_ext[I]*grad_u_ext[I]; */
1034  /* norm = sqrt(norm); */
1035  /* for(int I=0;I<nSpace;I++) */
1036  /* dir[I] = grad_u_ext[I]/norm; */
1037 
1038  /* ck.calculateGScale(G,dir,h_phi); */
1039  /* epsilon_redist = epsFact_redist*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter.data()[eN]); */
1040 
1041  /* // */
1042  /* //load the boundary values */
1043  /* // */
1044  /* bc_u_ext = isDOFBoundary_u.data()[ebNE_kb]*ebqe_bc_u_ext.data()[ebNE_kb]+(1-isDOFBoundary_u.data()[ebNE_kb])*u_ext; */
1045  /* // */
1046  /* //calculate the internal and external trace of the pde coefficients */
1047  /* // */
1048  /* evaluateCoefficients(epsilon_redist, */
1049  /* ebqe_phi_ls_ext.data()[ebNE_kb], */
1050  /* u_ext, */
1051  /* grad_u_ext, */
1052  /* m_ext, */
1053  /* dm_ext, */
1054  /* H_ext, */
1055  /* dH_ext, */
1056  /* r_ext); */
1057  /* evaluateCoefficients(epsilon_redist, */
1058  /* ebqe_phi_ls_ext.data()[ebNE_kb], */
1059  /* bc_u_ext, */
1060  /* bc_grad_u_ext, */
1061  /* bc_m_ext, */
1062  /* bc_dm_ext, */
1063  /* bc_H_ext, */
1064  /* bc_dH_ext, */
1065  /* bc_r_ext); */
1066  /* // */
1067  /* //calculate the numerical fluxes */
1068  /* // */
1069  /* //DoNothing for now */
1070  /* // */
1071  /* //calculate the flux jacobian */
1072  /* // */
1073  /* for (int j=0;j<nDOF_trial_element;j++) */
1074  /* { */
1075  /* //register int ebNE_kb_j = ebNE_kb*nDOF_trial_element+j; */
1076  /* //register int ebNE_kb_j_nSpace = ebNE_kb_j*nSpace; */
1077  /* //register int j_nSpace = j*nSpace; */
1078  /* register int ebN_local_kb_j=ebN_local_kb*nDOF_trial_element+j; */
1079 
1080  /* fluxJacobian_u_u[j]=ck.ExteriorNumericalAdvectiveFluxJacobian(dflux_u_u_ext,u_trial_trace_ref.data()[ebN_local_kb_j]); */
1081  /* }//j */
1082  /* // */
1083  /* //update the global Jacobian from the flux Jacobian */
1084  /* // */
1085  /* for (int i=0;i<nDOF_test_element;i++) */
1086  /* { */
1087  /* register int eN_i = eN*nDOF_test_element+i; */
1088  /* //register int ebNE_kb_i = ebNE_kb*nDOF_test_element+i; */
1089  /* for (int j=0;j<nDOF_trial_element;j++) */
1090  /* { */
1091  /* register int ebN_i_j = ebN*4*nDOF_test_X_trial_element + i*nDOF_trial_element + j; */
1092  /* //mwf debug */
1093  /* assert(fluxJacobian_u_u[j] == 0.0); */
1094  /* globalJacobian.data()[csrRowIndeces_u_u.data()[eN_i] + csrColumnOffsets_eb_u_u.data()[ebN_i_j]] += fluxJacobian_u_u[j]*u_test_dS[i]; */
1095  /* }//j */
1096  /* }//i */
1097  /* }//kb */
1098  /* }//ebNE */
1099  /* gf.useExact = useExact;//just to be safe */
1100  }//computeJacobian
1101 
1102  void calculateResidual_ellipticRedist(arguments_dict& args,
1103  bool useExact)
1104  {
1105  xt::pyarray<double>& mesh_trial_ref = args.array<double>("mesh_trial_ref");
1106  xt::pyarray<double>& mesh_grad_trial_ref = args.array<double>("mesh_grad_trial_ref");
1107  xt::pyarray<double>& mesh_dof = args.array<double>("mesh_dof");
1108  xt::pyarray<int>& mesh_l2g = args.array<int>("mesh_l2g");
1109  xt::pyarray<double>& x_ref = args.array<double>("x_ref");
1110  xt::pyarray<double>& dV_ref = args.array<double>("dV_ref");
1111  xt::pyarray<double>& u_trial_ref = args.array<double>("u_trial_ref");
1112  xt::pyarray<double>& u_grad_trial_ref = args.array<double>("u_grad_trial_ref");
1113  xt::pyarray<double>& u_test_ref = args.array<double>("u_test_ref");
1114  xt::pyarray<double>& u_grad_test_ref = args.array<double>("u_grad_test_ref");
1115  xt::pyarray<double>& mesh_trial_trace_ref = args.array<double>("mesh_trial_trace_ref");
1116  xt::pyarray<double>& mesh_grad_trial_trace_ref = args.array<double>("mesh_grad_trial_trace_ref");
1117  xt::pyarray<double>& dS_ref = args.array<double>("dS_ref");
1118  xt::pyarray<double>& u_trial_trace_ref = args.array<double>("u_trial_trace_ref");
1119  xt::pyarray<double>& u_grad_trial_trace_ref = args.array<double>("u_grad_trial_trace_ref");
1120  xt::pyarray<double>& u_test_trace_ref = args.array<double>("u_test_trace_ref");
1121  xt::pyarray<double>& u_grad_test_trace_ref = args.array<double>("u_grad_test_trace_ref");
1122  xt::pyarray<double>& normal_ref = args.array<double>("normal_ref");
1123  xt::pyarray<double>& boundaryJac_ref = args.array<double>("boundaryJac_ref");
1124  int nElements_global = args.scalar<int>("nElements_global");
1125  double useMetrics = args.scalar<double>("useMetrics");
1126  double alphaBDF = args.scalar<double>("alphaBDF");
1127  double epsFact_redist = args.scalar<double>("epsFact_redist");
1128  double backgroundDiffusionFactor = args.scalar<double>("backgroundDiffusionFactor");
1129  double weakDirichletFactor = args.scalar<double>("weakDirichletFactor");
1130  int freezeLevelSet = args.scalar<int>("freezeLevelSet");
1131  int useTimeIntegration = args.scalar<int>("useTimeIntegration");
1132  int lag_shockCapturing = args.scalar<int>("lag_shockCapturing");
1133  int lag_subgridError = args.scalar<int>("lag_subgridError");
1134  double shockCapturingDiffusion = args.scalar<double>("shockCapturingDiffusion");
1135  xt::pyarray<int>& u_l2g = args.array<int>("u_l2g");
1136  xt::pyarray<double>& elementDiameter = args.array<double>("elementDiameter");
1137  xt::pyarray<double>& nodeDiametersArray = args.array<double>("nodeDiametersArray");
1138  xt::pyarray<double>& u_dof = args.array<double>("u_dof");
1139  xt::pyarray<double>& phi_dof = args.array<double>("phi_dof");
1140  xt::pyarray<double>& phi_ls = args.array<double>("phi_ls");
1141  xt::pyarray<double>& q_m = args.array<double>("q_m");
1142  xt::pyarray<double>& q_u = args.array<double>("q_u");
1143  xt::pyarray<double>& q_n = args.array<double>("q_n");
1144  xt::pyarray<double>& q_dH = args.array<double>("q_dH");
1145  xt::pyarray<double>& u_weak_internal_bc_dofs = args.array<double>("u_weak_internal_bc_dofs");
1146  xt::pyarray<double>& q_m_betaBDF = args.array<double>("q_m_betaBDF");
1147  xt::pyarray<double>& q_dH_last = args.array<double>("q_dH_last");
1148  xt::pyarray<double>& q_cfl = args.array<double>("q_cfl");
1149  xt::pyarray<double>& q_numDiff_u = args.array<double>("q_numDiff_u");
1150  xt::pyarray<double>& q_numDiff_u_last = args.array<double>("q_numDiff_u_last");
1151  xt::pyarray<int>& weakDirichletConditionFlags = args.array<int>("weakDirichletConditionFlags");
1152  int offset_u = args.scalar<int>("offset_u");
1153  int stride_u = args.scalar<int>("stride_u");
1154  xt::pyarray<double>& globalResidual = args.array<double>("globalResidual");
1155  int nExteriorElementBoundaries_global = args.scalar<int>("nExteriorElementBoundaries_global");
1156  xt::pyarray<int>& exteriorElementBoundariesArray = args.array<int>("exteriorElementBoundariesArray");
1157  xt::pyarray<int>& elementBoundaryElementsArray = args.array<int>("elementBoundaryElementsArray");
1158  xt::pyarray<int>& elementBoundaryLocalElementBoundariesArray = args.array<int>("elementBoundaryLocalElementBoundariesArray");
1159  xt::pyarray<double>& ebqe_phi_ls_ext = args.array<double>("ebqe_phi_ls_ext");
1160  xt::pyarray<int>& isDOFBoundary_u = args.array<int>("isDOFBoundary_u");
1161  xt::pyarray<double>& ebqe_bc_u_ext = args.array<double>("ebqe_bc_u_ext");
1162  xt::pyarray<double>& ebqe_u = args.array<double>("ebqe_u");
1163  xt::pyarray<double>& ebqe_n = args.array<double>("ebqe_n");
1164  int ELLIPTIC_REDISTANCING = args.scalar<int>("ELLIPTIC_REDISTANCING");
1165  double backgroundDissipationEllipticRedist = args.scalar<double>("backgroundDissipationEllipticRedist");
1166  xt::pyarray<double>& lumped_qx = args.array<double>("lumped_qx");
1167  xt::pyarray<double>& lumped_qy = args.array<double>("lumped_qy");
1168  xt::pyarray<double>& lumped_qz = args.array<double>("lumped_qz");
1169  double alpha = args.scalar<double>("alpha");
1170  gf.useExact=useExact;
1171  //
1172  //loop over elements to compute volume integrals and load them into element and global residual
1173  //
1174  //eN is the element index
1175  //eN_k is the quadrature point index for a scalar
1176  //eN_k_nSpace is the quadrature point index for a vector
1177  //eN_i is the element test function index
1178  //eN_j is the element trial function index
1179  //eN_k_j is the quadrature point index for a trial function
1180  //eN_k_i is the quadrature point index for a trial function
1181  for(int eN=0;eN<nElements_global;eN++)
1182  {
1183  //declare local storage for element residual and initialize
1184  register double elementResidual_u[nDOF_test_element],element_phi[nDOF_trial_element];
1185  double epsilon_redist,h_phi, norm;
1186  for (int i=0;i<nDOF_test_element;i++)
1187  {
1188  register int eN_i=eN*nDOF_trial_element+i;
1189  elementResidual_u[i]=0.0;
1190  element_phi[i] = phi_dof.data()[u_l2g.data()[eN_i]];
1191  }//i
1192  double element_nodes[nDOF_mesh_trial_element*3];
1193  for (int i=0;i<nDOF_mesh_trial_element;i++)
1194  {
1195  register int eN_i=eN*nDOF_mesh_trial_element+i;
1196  for(int I=0;I<3;I++)
1197  element_nodes[i*3 + I] = mesh_dof.data()[mesh_l2g.data()[eN_i]*3 + I];
1198  }//i
1199  gf.calculate(element_phi, element_nodes, x_ref.data(),false);
1200  //loop over quadrature points and compute integrands
1201  for (int k=0;k<nQuadraturePoints_element;k++)
1202  {
1203  gf.set_quad(k);
1204  //compute indeces and declare local storage
1205  register int eN_k = eN*nQuadraturePoints_element+k,
1206  eN_k_nSpace = eN_k*nSpace,
1207  eN_nDOF_trial_element = eN*nDOF_trial_element;
1208  register double
1209  coeff, delta,
1210  qx, qy, qz, normalReconstruction[nSpace],
1211  u=0,grad_u[nSpace],
1212  m=0.0,
1213  jac[nSpace*nSpace], jacDet, jacInv[nSpace*nSpace],
1214  u_grad_trial[nDOF_trial_element*nSpace],
1215  u_test_dV[nDOF_trial_element], u_grad_test_dV[nDOF_test_element*nSpace],
1216  dV,x,y,z,G[nSpace*nSpace],G_dd_G,tr_G;
1217  ck.calculateMapping_element(eN,
1218  k,
1219  mesh_dof.data(),
1220  mesh_l2g.data(),
1221  mesh_trial_ref.data(),
1222  mesh_grad_trial_ref.data(),
1223  jac,
1224  jacDet,
1225  jacInv,
1226  x,y,z);
1227  ck.calculateH_element(eN,
1228  k,
1229  nodeDiametersArray.data(),
1230  mesh_l2g.data(),
1231  mesh_trial_ref.data(),
1232  h_phi);
1233  //get the physical integration weight
1234  dV = fabs(jacDet)*dV_ref.data()[k];
1235  ck.calculateG(jacInv,G,G_dd_G,tr_G);
1236  //get the trial function gradients
1237  ck.gradTrialFromRef(&u_grad_trial_ref.data()[k*nDOF_trial_element*nSpace],
1238  jacInv,
1239  u_grad_trial);
1240  //get the solution
1241  ck.valFromDOF(u_dof.data(),
1242  &u_l2g.data()[eN_nDOF_trial_element],
1243  &u_trial_ref.data()[k*nDOF_trial_element],
1244  u);
1245  //get the solution gradients
1246  ck.gradFromDOF(u_dof.data(),
1247  &u_l2g.data()[eN_nDOF_trial_element],
1248  u_grad_trial,
1249  grad_u);
1250  if (ELLIPTIC_REDISTANCING > 1)
1251  { // use linear elliptic re-distancing via C0 normal reconstruction
1252  ck.valFromDOF(lumped_qx.data(),
1253  &u_l2g.data()[eN_nDOF_trial_element],
1254  &u_trial_ref.data()[k*nDOF_trial_element],
1255  qx);
1256  ck.valFromDOF(lumped_qy.data(),
1257  &u_l2g.data()[eN_nDOF_trial_element],
1258  &u_trial_ref.data()[k*nDOF_trial_element],
1259  qy);
1260  ck.valFromDOF(lumped_qz.data(),
1261  &u_l2g.data()[eN_nDOF_trial_element],
1262  &u_trial_ref.data()[k*nDOF_trial_element],
1263  qz);
1264  normalReconstruction[0] = qx;
1265  normalReconstruction[1] = qy;
1266  if (nSpace == 3)
1267  normalReconstruction[2] = qz;
1268  }
1269  //precalculate test function products with integration weights
1270  for (int j=0;j<nDOF_trial_element;j++)
1271  {
1272  u_test_dV[j] = u_test_ref.data()[k*nDOF_trial_element+j]*dV;
1273  for (int I=0;I<nSpace;I++)
1274  u_grad_test_dV[j*nSpace+I] = u_grad_trial[j*nSpace+I]*dV;
1275  }
1276  // MOVING MESH. Omit for now //
1277  // COMPUTE NORM OF GRAD(u) //
1278  double norm_grad_u = 0.;
1279  for (int I=0;I<nSpace;I++)
1280  norm_grad_u += grad_u[I]*grad_u[I];
1281  norm_grad_u = std::sqrt(norm_grad_u) + 1.0E-10;
1282 
1283  // SAVE MASS AND SOLUTION FOR OTHER MODELS //
1284  q_m.data()[eN_k] = u; //m=u
1285  q_u.data()[eN_k] = u;
1286  for (int I=0;I<nSpace;I++)
1287  q_n.data()[eN_k_nSpace+I] = grad_u[I]/norm_grad_u;
1288 
1289  // COMPUTE COEFFICIENTS //
1290  if (SINGLE_POTENTIAL == 1)
1291  coeff = 1.0-1.0/norm_grad_u; //single potential
1292  else // double potential
1293  coeff = 1.0+2*std::pow(norm_grad_u,2)-3*norm_grad_u;
1294 
1295  // COMPUTE DELTA FUNCTION //
1296  epsilon_redist = epsFact_redist*(useMetrics*h_phi
1297  +(1.0-useMetrics)*elementDiameter.data()[eN]);
1298  delta = gf.D(epsilon_redist,phi_ls.data()[eN_k]);
1299 
1300  // COMPUTE STRONG RESIDUAL //
1301  double Si = -1.0+2.0*gf.H(epsilon_redist,phi_ls.data()[eN_k]);
1302  double residualEikonal = Si*(norm_grad_u-1.0);
1303  double backgroundDissipation = backgroundDissipationEllipticRedist*elementDiameter.data()[eN];
1304 
1305  // UPDATE ELEMENT RESIDUAL //
1306  for(int i=0;i<nDOF_test_element;i++)
1307  {
1308  register int i_nSpace = i*nSpace;
1309  // global i-th index
1310  int gi = offset_u+stride_u*u_l2g.data()[eN*nDOF_test_element+i];
1311 
1312  if (ELLIPTIC_REDISTANCING > 1) // (NON)LINEAR VIA C0 NORMAL RECONSTRUCTION
1313  {
1314  elementResidual_u[i] +=
1315  residualEikonal*u_test_dV[i]
1316  +ck.NumericalDiffusion(1.0+backgroundDissipation,
1317  grad_u,
1318  &u_grad_test_dV[i_nSpace])
1319  -ck.NumericalDiffusion(1.0,
1320  normalReconstruction,
1321  &u_grad_test_dV[i_nSpace])
1322  +alpha*(u_dof.data()[gi]-phi_dof.data()[gi])*delta*u_test_dV[i]; // BCs
1323  }
1324  else // =1. Nonlinear via single or double pot.
1325  {
1326  elementResidual_u[i] +=
1327  residualEikonal*u_test_dV[i]
1328  +ck.NumericalDiffusion(coeff+backgroundDissipation,
1329  grad_u,
1330  &u_grad_test_dV[i_nSpace])
1331  + alpha*(u_dof.data()[gi]-phi_dof.data()[gi])*delta*u_test_dV[i]; // BCs
1332  }
1333  }//i
1334  }//k
1335  //
1336  //load element into global residual and save element residual
1337  //
1338  for(int i=0;i<nDOF_test_element;i++)
1339  {
1340  register int eN_i=eN*nDOF_test_element+i;
1341  globalResidual.data()[offset_u+stride_u*u_l2g.data()[eN_i]]+=elementResidual_u[i];
1342  }//i
1343  }//elements
1344  //
1345  //loop over exterior element boundaries to save soln at quad points
1346  //
1347  for (int ebNE = 0; ebNE < nExteriorElementBoundaries_global; ebNE++)
1348  {
1349  register int ebN = exteriorElementBoundariesArray.data()[ebNE],
1350  eN = elementBoundaryElementsArray.data()[ebN*2+0],
1351  ebN_local = elementBoundaryLocalElementBoundariesArray.data()[ebN*2+0],
1352  eN_nDOF_trial_element = eN*nDOF_trial_element;
1353  for (int kb=0;kb<nQuadraturePoints_elementBoundary;kb++)
1354  {
1355  register int ebNE_kb = ebNE*nQuadraturePoints_elementBoundary+kb,
1356  ebNE_kb_nSpace = ebNE_kb*nSpace,
1357  ebN_local_kb = ebN_local*nQuadraturePoints_elementBoundary+kb,
1358  ebN_local_kb_nSpace = ebN_local_kb*nSpace;
1359  register double
1360  u_ext=0.0,
1361  grad_u_ext[nSpace],
1362  jac_ext[nSpace*nSpace],jacDet_ext,jacInv_ext[nSpace*nSpace],
1363  boundaryJac[nSpace*(nSpace-1)],
1364  metricTensor[(nSpace-1)*(nSpace-1)],metricTensorDetSqrt,
1365  u_grad_trial_trace[nDOF_trial_element*nSpace],
1366  normal[nSpace],x_ext,y_ext,z_ext,
1367  dir[nSpace],norm;
1368  ck.calculateMapping_elementBoundary(eN,
1369  ebN_local,
1370  kb,
1371  ebN_local_kb,
1372  mesh_dof.data(),
1373  mesh_l2g.data(),
1374  mesh_trial_trace_ref.data(),
1375  mesh_grad_trial_trace_ref.data(),
1376  boundaryJac_ref.data(),
1377  jac_ext,
1378  jacDet_ext,
1379  jacInv_ext,
1380  boundaryJac,
1381  metricTensor,
1382  metricTensorDetSqrt,
1383  normal_ref.data(),
1384  normal,
1385  x_ext,y_ext,z_ext);
1386  //compute shape and solution information
1387  //shape
1388  ck.gradTrialFromRef(&u_grad_trial_trace_ref.data()[ebN_local_kb_nSpace*nDOF_trial_element],
1389  jacInv_ext,
1390  u_grad_trial_trace);
1391  //solution and gradients
1392  ck.valFromDOF(u_dof.data(),
1393  &u_l2g.data()[eN_nDOF_trial_element],
1394  &u_trial_trace_ref.data()[ebN_local_kb*nDOF_test_element],
1395  u_ext);
1396  ck.gradFromDOF(u_dof.data(),
1397  &u_l2g.data()[eN_nDOF_trial_element],
1398  u_grad_trial_trace,
1399  grad_u_ext);
1400  norm = 0;
1401  for (int I=0;I<nSpace;I++)
1402  norm += grad_u_ext[I]*grad_u_ext[I];
1403  norm = sqrt(norm) + 1.0E-10;
1404  for (int I=0;I<nSpace;I++)
1405  dir[I] = grad_u_ext[I]/norm;
1406 
1407  //save for other models
1408  ebqe_u.data()[ebNE_kb] = u_ext;
1409 
1410  for (int I=0;I<nSpace;I++)
1411  ebqe_n.data()[ebNE_kb_nSpace+I] = dir[I];
1412  }//kb
1413  }//ebNE
1414  }
1415 
1416  void calculateJacobian_ellipticRedist(arguments_dict& args,
1417  bool useExact)
1418  {
1419  xt::pyarray<double>& mesh_trial_ref = args.array<double>("mesh_trial_ref");
1420  xt::pyarray<double>& mesh_grad_trial_ref = args.array<double>("mesh_grad_trial_ref");
1421  xt::pyarray<double>& mesh_dof = args.array<double>("mesh_dof");
1422  xt::pyarray<int>& mesh_l2g = args.array<int>("mesh_l2g");
1423  xt::pyarray<double>& x_ref = args.array<double>("x_ref");
1424  xt::pyarray<double>& dV_ref = args.array<double>("dV_ref");
1425  xt::pyarray<double>& u_trial_ref = args.array<double>("u_trial_ref");
1426  xt::pyarray<double>& u_grad_trial_ref = args.array<double>("u_grad_trial_ref");
1427  xt::pyarray<double>& u_test_ref = args.array<double>("u_test_ref");
1428  xt::pyarray<double>& u_grad_test_ref = args.array<double>("u_grad_test_ref");
1429  xt::pyarray<double>& mesh_trial_trace_ref = args.array<double>("mesh_trial_trace_ref");
1430  xt::pyarray<double>& mesh_grad_trial_trace_ref = args.array<double>("mesh_grad_trial_trace_ref");
1431  xt::pyarray<double>& dS_ref = args.array<double>("dS_ref");
1432  xt::pyarray<double>& u_trial_trace_ref = args.array<double>("u_trial_trace_ref");
1433  xt::pyarray<double>& u_grad_trial_trace_ref = args.array<double>("u_grad_trial_trace_ref");
1434  xt::pyarray<double>& u_test_trace_ref = args.array<double>("u_test_trace_ref");
1435  xt::pyarray<double>& u_grad_test_trace_ref = args.array<double>("u_grad_test_trace_ref");
1436  xt::pyarray<double>& normal_ref = args.array<double>("normal_ref");
1437  xt::pyarray<double>& boundaryJac_ref = args.array<double>("boundaryJac_ref");
1438  int nElements_global = args.scalar<int>("nElements_global");
1439  double useMetrics = args.scalar<double>("useMetrics");
1440  double alphaBDF = args.scalar<double>("alphaBDF");
1441  double epsFact_redist = args.scalar<double>("epsFact_redist");
1442  double backgroundDiffusionFactor = args.scalar<double>("backgroundDiffusionFactor");
1443  double weakDirichletFactor = args.scalar<double>("weakDirichletFactor");
1444  int freezeLevelSet = args.scalar<int>("freezeLevelSet");
1445  int useTimeIntegration = args.scalar<int>("useTimeIntegration");
1446  int lag_shockCapturing = args.scalar<int>("lag_shockCapturing");
1447  int lag_subgridError = args.scalar<int>("lag_subgridError");
1448  double shockCapturingDiffusion = args.scalar<double>("shockCapturingDiffusion");
1449  xt::pyarray<int>& u_l2g = args.array<int>("u_l2g");
1450  xt::pyarray<double>& elementDiameter = args.array<double>("elementDiameter");
1451  xt::pyarray<double>& nodeDiametersArray = args.array<double>("nodeDiametersArray");
1452  xt::pyarray<double>& u_dof = args.array<double>("u_dof");
1453  xt::pyarray<double>& phi_dof = args.array<double>("phi_dof");
1454  xt::pyarray<double>& u_weak_internal_bc_dofs = args.array<double>("u_weak_internal_bc_dofs");
1455  xt::pyarray<double>& phi_ls = args.array<double>("phi_ls");
1456  xt::pyarray<double>& q_m_betaBDF = args.array<double>("q_m_betaBDF");
1457  xt::pyarray<double>& q_dH_last = args.array<double>("q_dH_last");
1458  xt::pyarray<double>& q_cfl = args.array<double>("q_cfl");
1459  xt::pyarray<double>& q_numDiff_u = args.array<double>("q_numDiff_u");
1460  xt::pyarray<double>& q_numDiff_u_last = args.array<double>("q_numDiff_u_last");
1461  xt::pyarray<int>& weakDirichletConditionFlags = args.array<int>("weakDirichletConditionFlags");
1462  xt::pyarray<int>& csrRowIndeces_u_u = args.array<int>("csrRowIndeces_u_u");
1463  xt::pyarray<int>& csrColumnOffsets_u_u = args.array<int>("csrColumnOffsets_u_u");
1464  xt::pyarray<double>& globalJacobian = args.array<double>("globalJacobian");
1465  int nExteriorElementBoundaries_global = args.scalar<int>("nExteriorElementBoundaries_global");
1466  xt::pyarray<int>& exteriorElementBoundariesArray = args.array<int>("exteriorElementBoundariesArray");
1467  xt::pyarray<int>& elementBoundaryElementsArray = args.array<int>("elementBoundaryElementsArray");
1468  xt::pyarray<int>& elementBoundaryLocalElementBoundariesArray = args.array<int>("elementBoundaryLocalElementBoundariesArray");
1469  xt::pyarray<double>& ebqe_phi_ls_ext = args.array<double>("ebqe_phi_ls_ext");
1470  xt::pyarray<int>& isDOFBoundary_u = args.array<int>("isDOFBoundary_u");
1471  xt::pyarray<double>& ebqe_bc_u_ext = args.array<double>("ebqe_bc_u_ext");
1472  xt::pyarray<int>& csrColumnOffsets_eb_u_u = args.array<int>("csrColumnOffsets_eb_u_u");
1473  int ELLIPTIC_REDISTANCING = args.scalar<int>("ELLIPTIC_REDISTANCING");
1474  double backgroundDissipationEllipticRedist = args.scalar<double>("backgroundDissipationEllipticRedist");
1475  double alpha = args.scalar<double>("alpha");
1476  gf.useExact=useExact;
1477  //
1478  //loop over elements
1479  //
1480  for(int eN=0;eN<nElements_global;eN++)
1481  {
1482  register double elementJacobian_u_u[nDOF_test_element][nDOF_trial_element],element_phi[nDOF_trial_element];
1483  double epsilon_redist,h_phi, norm;
1484  for (int i=0;i<nDOF_test_element;i++)
1485  {
1486  register int eN_i=eN*nDOF_trial_element+i;
1487  element_phi[i] = phi_dof.data()[u_l2g.data()[eN_i]];
1488  for (int j=0;j<nDOF_trial_element;j++)
1489  {
1490  elementJacobian_u_u[i][j]=0.0;
1491  }
1492  }
1493  double element_nodes[nDOF_mesh_trial_element*3];
1494  for (int i=0;i<nDOF_mesh_trial_element;i++)
1495  {
1496  register int eN_i=eN*nDOF_mesh_trial_element+i;
1497  for(int I=0;I<3;I++)
1498  element_nodes[i*3 + I] = mesh_dof.data()[mesh_l2g.data()[eN_i]*3 + I];
1499  }//i
1500  gf.calculate(element_phi, element_nodes, x_ref.data(),false);
1501  for (int k=0;k<nQuadraturePoints_element;k++)
1502  {
1503  gf.set_quad(k);
1504  int eN_k = eN*nQuadraturePoints_element+k, //index to a scalar at a quadrature point
1505  eN_k_nSpace = eN_k*nSpace,
1506  eN_nDOF_trial_element = eN*nDOF_trial_element; //index to a vector at a quadrature point
1507 
1508  //declare local storage
1509  register double
1510  coeff1, coeff2, delta,
1511  grad_u[nSpace],
1512  jac[nSpace*nSpace], jacDet, jacInv[nSpace*nSpace],
1513  u_grad_trial[nDOF_trial_element*nSpace],
1514  dV, u_test_dV[nDOF_test_element], u_grad_test_dV[nDOF_test_element*nSpace],
1515  x,y,z,G[nSpace*nSpace],G_dd_G,tr_G;
1516  //
1517  //calculate solution and gradients at quadrature points
1518  //
1519  ck.calculateMapping_element(eN,
1520  k,
1521  mesh_dof.data(),
1522  mesh_l2g.data(),
1523  mesh_trial_ref.data(),
1524  mesh_grad_trial_ref.data(),
1525  jac,
1526  jacDet,
1527  jacInv,
1528  x,y,z);
1529  ck.calculateH_element(eN,
1530  k,
1531  nodeDiametersArray.data(),
1532  mesh_l2g.data(),
1533  mesh_trial_ref.data(),
1534  h_phi);
1535  //get the physical integration weight
1536  dV = fabs(jacDet)*dV_ref.data()[k];
1537  ck.calculateG(jacInv,G,G_dd_G,tr_G);
1538  //get the trial function gradients
1539  ck.gradTrialFromRef(&u_grad_trial_ref.data()[k*nDOF_trial_element*nSpace],
1540  jacInv,
1541  u_grad_trial);
1542  //get the solution gradients
1543  ck.gradFromDOF(u_dof.data(),
1544  &u_l2g.data()[eN_nDOF_trial_element],
1545  u_grad_trial,
1546  grad_u);
1547  //precalculate test function products with integration weights
1548  for (int j=0;j<nDOF_trial_element;j++)
1549  {
1550  u_test_dV[j] = u_test_ref.data()[k*nDOF_trial_element+j]*dV;
1551  for (int I=0;I<nSpace;I++)
1552  u_grad_test_dV[j*nSpace+I] = u_grad_trial[j*nSpace+I]*dV;
1553  }
1554  // MOVING MESH. Omit for now //
1555  // COMPUTE NORM OF GRAD(u) //
1556  double norm_grad_u = 0;
1557  for(int I=0;I<nSpace;I++)
1558  norm_grad_u += grad_u[I]*grad_u[I];
1559  norm_grad_u = std::sqrt(norm_grad_u) + 1.0E-10;
1560 
1561  // COMPUTE COEFFICIENTS //
1562  if (SINGLE_POTENTIAL==1)
1563  {
1564  coeff1 = 0.; //-1./norm_grad_u;
1565  coeff2 = 1./std::pow(norm_grad_u,3);
1566  }
1567  else
1568  {
1569  coeff1 = fmax(1.0E-10, 2*std::pow(norm_grad_u,2)-3*norm_grad_u);
1570  coeff2 = fmax(1.0E-10, 4.-3./norm_grad_u);
1571  }
1572 
1573  // COMPUTE DELTA FUNCTION //
1574  epsilon_redist = epsFact_redist*(useMetrics*h_phi
1575  +(1.0-useMetrics)*elementDiameter.data()[eN]);
1576  delta = gf.D(epsilon_redist,phi_ls.data()[eN_k]);
1577 
1578  // COMPUTE STRONG Jacobian //
1579  double Si = -1.0+2.0*gf.H(epsilon_redist,phi_ls.data()[eN_k]);
1580  double dH[nSpace];
1581  for (int I=0; I<nSpace;I++)
1582  dH[I] = Si*grad_u[I]/norm_grad_u;
1583  double backgroundDissipation = backgroundDissipationEllipticRedist*elementDiameter.data()[eN];
1584 
1585  // LOOP IN I-DOFs //
1586  for(int i=0;i<nDOF_test_element;i++)
1587  {
1588  int i_nSpace = i*nSpace;
1589  for(int j=0;j<nDOF_trial_element;j++)
1590  {
1591  int j_nSpace = j*nSpace;
1592  elementJacobian_u_u[i][j] +=
1593  ck.HamiltonianJacobian_weak(dH,&u_grad_trial[j_nSpace],u_test_dV[i])
1594  +ck.NumericalDiffusionJacobian(1.0+backgroundDissipation,
1595  &u_grad_trial[j_nSpace],
1596  &u_grad_test_dV[i_nSpace])
1597  + (ELLIPTIC_REDISTANCING == 1 ? 1. : 0.)*
1598  ( ck.NumericalDiffusionJacobian(coeff1,
1599  &u_grad_trial[j_nSpace],
1600  &u_grad_test_dV[i_nSpace])
1601  + coeff2*dV*
1602  ck.NumericalDiffusion(1.0,grad_u,&u_grad_trial[i_nSpace])*
1603  ck.NumericalDiffusion(1.0,grad_u,&u_grad_trial[j_nSpace]) )
1604  + (i == j ? alpha*delta*u_test_dV[i] : 0.); //lumped
1605  }//j
1606  }//i
1607  }//k
1608  //
1609  //load into element Jacobian into global Jacobian
1610  //
1611  for (int i=0;i<nDOF_test_element;i++)
1612  {
1613  int eN_i = eN*nDOF_test_element+i;
1614  for (int j=0;j<nDOF_trial_element;j++)
1615  {
1616  int eN_i_j = eN_i*nDOF_trial_element+j;
1617  globalJacobian.data()[csrRowIndeces_u_u.data()[eN_i]
1618  + csrColumnOffsets_u_u.data()[eN_i_j]] += elementJacobian_u_u[i][j];
1619  }//j
1620  }//i
1621  }//elements
1622  }//computeJacobian
1623 
1624  void normalReconstruction(arguments_dict& args)
1625  {
1626  xt::pyarray<double>& mesh_trial_ref = args.array<double>("mesh_trial_ref");
1627  xt::pyarray<double>& mesh_grad_trial_ref = args.array<double>("mesh_grad_trial_ref");
1628  xt::pyarray<double>& mesh_dof = args.array<double>("mesh_dof");
1629  xt::pyarray<int>& mesh_l2g = args.array<int>("mesh_l2g");
1630  xt::pyarray<double>& dV_ref = args.array<double>("dV_ref");
1631  xt::pyarray<double>& u_trial_ref = args.array<double>("u_trial_ref");
1632  xt::pyarray<double>& u_grad_trial_ref = args.array<double>("u_grad_trial_ref");
1633  xt::pyarray<double>& u_test_ref = args.array<double>("u_test_ref");
1634  int nElements_global = args.scalar<int>("nElements_global");
1635  xt::pyarray<int>& u_l2g = args.array<int>("u_l2g");
1636  xt::pyarray<double>& elementDiameter = args.array<double>("elementDiameter");
1637  xt::pyarray<double>& u_dof = args.array<double>("u_dof");
1638  int offset_u = args.scalar<int>("offset_u");
1639  int stride_u = args.scalar<int>("stride_u");
1640  int numDOFs = args.scalar<int>("numDOFs");
1641  xt::pyarray<double>& lumped_qx = args.array<double>("lumped_qx");
1642  xt::pyarray<double>& lumped_qy = args.array<double>("lumped_qy");
1643  xt::pyarray<double>& lumped_qz = args.array<double>("lumped_qz");
1644  weighted_lumped_mass_matrix.resize(numDOFs,0.0);
1645  for (int i=0; i<numDOFs; i++)
1646  {
1647  // output vectors
1648  lumped_qx.data()[i]=0.;
1649  lumped_qy.data()[i]=0.;
1650  lumped_qz.data()[i]=0.;
1651  // auxiliary vectors
1652  weighted_lumped_mass_matrix[i]=0.;
1653  }
1654  for(int eN=0;eN<nElements_global;eN++)
1655  {
1656  //declare local storage for local contributions and initialize
1657  register double
1658  element_weighted_lumped_mass_matrix[nDOF_test_element],
1659  element_rhsx_normal_reconstruction[nDOF_test_element],
1660  element_rhsy_normal_reconstruction[nDOF_test_element],
1661  element_rhsz_normal_reconstruction[nDOF_test_element];
1662  for (int i=0;i<nDOF_test_element;i++)
1663  {
1664  element_weighted_lumped_mass_matrix[i]=0.0;
1665  element_rhsx_normal_reconstruction[i]=0.0;
1666  element_rhsy_normal_reconstruction[i]=0.0;
1667  element_rhsz_normal_reconstruction[i]=0.0;
1668  }
1669  //loop over quadrature points and compute integrands
1670  for (int k=0;k<nQuadraturePoints_element;k++)
1671  {
1672  //compute indeces and declare local storage
1673  register int eN_k = eN*nQuadraturePoints_element+k,
1674  eN_k_nSpace = eN_k*nSpace,
1675  eN_nDOF_trial_element = eN*nDOF_trial_element;
1676  register double
1677  //for mass matrix contributions
1678  grad_u[nSpace],
1679  u_grad_trial[nDOF_trial_element*nSpace],
1680  u_test_dV[nDOF_trial_element],
1681  //for general use
1682  jac[nSpace*nSpace], jacDet, jacInv[nSpace*nSpace],
1683  dV,x,y,z;
1684  //get the physical integration weight
1685  ck.calculateMapping_element(eN,
1686  k,
1687  mesh_dof.data(),
1688  mesh_l2g.data(),
1689  mesh_trial_ref.data(),
1690  mesh_grad_trial_ref.data(),
1691  jac,
1692  jacDet,
1693  jacInv,
1694  x,y,z);
1695  dV = fabs(jacDet)*dV_ref.data()[k];
1696  ck.gradTrialFromRef(&u_grad_trial_ref.data()[k*nDOF_trial_element*nSpace],
1697  jacInv,
1698  u_grad_trial);
1699  ck.gradFromDOF(u_dof.data(),
1700  &u_l2g.data()[eN_nDOF_trial_element],u_grad_trial,
1701  grad_u);
1702  //precalculate test function products with integration weights for mass matrix terms
1703  for (int j=0;j<nDOF_trial_element;j++)
1704  u_test_dV[j] = u_test_ref.data()[k*nDOF_trial_element+j]*dV;
1705 
1706  double rhsx = grad_u[0];
1707  double rhsy = grad_u[1];
1708  double rhsz = 0;
1709  if (nSpace==3)
1710  rhsz = grad_u[2];
1711 
1712  double norm_grad_u = 0;
1713  for (int I=0;I<nSpace; I++)
1714  norm_grad_u += grad_u[I]*grad_u[I];
1715  norm_grad_u = std::sqrt(norm_grad_u) + 1.0E-10;
1716 
1717  for(int i=0;i<nDOF_test_element;i++)
1718  {
1719  element_weighted_lumped_mass_matrix[i] += norm_grad_u*u_test_dV[i];
1720  element_rhsx_normal_reconstruction[i] += rhsx*u_test_dV[i];
1721  element_rhsy_normal_reconstruction[i] += rhsy*u_test_dV[i];
1722  element_rhsz_normal_reconstruction[i] += rhsz*u_test_dV[i];
1723  }
1724  } //k
1725  // DISTRIBUTE //
1726  for(int i=0;i<nDOF_test_element;i++)
1727  {
1728  int eN_i=eN*nDOF_test_element+i;
1729  int gi = offset_u+stride_u*u_l2g.data()[eN_i]; //global i-th index
1730 
1731  weighted_lumped_mass_matrix[gi] += element_weighted_lumped_mass_matrix[i];
1732  lumped_qx.data()[gi] += element_rhsx_normal_reconstruction[i];
1733  lumped_qy.data()[gi] += element_rhsy_normal_reconstruction[i];
1734  lumped_qz.data()[gi] += element_rhsz_normal_reconstruction[i];
1735  }//i
1736  }//elements
1737  // COMPUTE LUMPED L2 PROJECTION
1738  for (int i=0; i<numDOFs; i++)
1739  {
1740  // normal reconstruction
1741  double weighted_mi = weighted_lumped_mass_matrix[i];
1742  lumped_qx.data()[i] /= weighted_mi;
1743  lumped_qy.data()[i] /= weighted_mi;
1744  lumped_qz.data()[i] /= weighted_mi;
1745  }
1746  }
1747 
1748  std::tuple<double, double, double> calculateMetricsAtEOS(arguments_dict& args)
1749  {
1750  xt::pyarray<double>& mesh_trial_ref = args.array<double>("mesh_trial_ref");
1751  xt::pyarray<double>& mesh_grad_trial_ref = args.array<double>("mesh_grad_trial_ref");
1752  xt::pyarray<double>& mesh_dof = args.array<double>("mesh_dof");
1753  xt::pyarray<int>& mesh_l2g = args.array<int>("mesh_l2g");
1754  xt::pyarray<double>& dV_ref = args.array<double>("dV_ref");
1755  xt::pyarray<double>& u_trial_ref = args.array<double>("u_trial_ref");
1756  xt::pyarray<double>& u_grad_trial_ref = args.array<double>("u_grad_trial_ref");
1757  xt::pyarray<double>& u_test_ref = args.array<double>("u_test_ref");
1758  int nElements_global = args.scalar<int>("nElements_global");
1759  xt::pyarray<int>& u_l2g = args.array<int>("u_l2g");
1760  xt::pyarray<double>& elementDiameter = args.array<double>("elementDiameter");
1761  double degree_polynomial = args.scalar<double>("degree_polynomial");
1762  double epsFact_redist = args.scalar<double>("epsFact_redist");
1763  xt::pyarray<double>& u_dof = args.array<double>("u_dof");
1764  xt::pyarray<double>& u_exact = args.array<double>("u_exact");
1765  int offset_u = args.scalar<int>("offset_u");
1766  int stride_u = args.scalar<int>("stride_u)");
1767  double global_V = 0.;
1768  double global_V0 = 0.;
1769  double global_I_err = 0.0;
1770  double global_V_err = 0.0;
1771  double global_D_err = 0.0;
1773  // ** LOOP IN CELLS //
1775  for(int eN=0;eN<nElements_global;eN++)
1776  {
1777  //declare local storage for local contributions and initialize
1778  register double
1779  elementResidual_u[nDOF_test_element];
1780  double
1781  cell_mass_error = 0., cell_mass_exact = 0.,
1782  cell_I_err = 0.,
1783  cell_V = 0., cell_V0 = 0.,
1784  cell_D_err = 0.;
1785 
1786  //loop over quadrature points and compute integrands
1787  for (int k=0;k<nQuadraturePoints_element;k++)
1788  {
1789  //compute indeces and declare local storage
1790  register int eN_k = eN*nQuadraturePoints_element+k,
1791  eN_k_nSpace = eN_k*nSpace,
1792  eN_nDOF_trial_element = eN*nDOF_trial_element;
1793  register double
1794  u, uh,
1795  u_grad_trial[nDOF_trial_element*nSpace],
1796  grad_uh[nSpace],
1797  //for general use
1798  jac[nSpace*nSpace], jacDet, jacInv[nSpace*nSpace],
1799  dV,x,y,z;
1800  //get the physical integration weight
1801  ck.calculateMapping_element(eN,
1802  k,
1803  mesh_dof.data(),
1804  mesh_l2g.data(),
1805  mesh_trial_ref.data(),
1806  mesh_grad_trial_ref.data(),
1807  jac,
1808  jacDet,
1809  jacInv,
1810  x,y,z);
1811  dV = fabs(jacDet)*dV_ref.data()[k];
1812  // get functions at quad points
1813  ck.valFromDOF(u_dof.data(),
1814  &u_l2g.data()[eN_nDOF_trial_element],&u_trial_ref.data()[k*nDOF_trial_element],
1815  uh);
1816  u = u_exact.data()[eN_k];
1817  // get gradients
1818  ck.gradTrialFromRef(&u_grad_trial_ref.data()[k*nDOF_trial_element*nSpace],
1819  jacInv,
1820  u_grad_trial);
1821  ck.gradFromDOF(u_dof.data(),&u_l2g.data()[eN_nDOF_trial_element],u_grad_trial,grad_uh);
1822 
1823  double epsHeaviside = epsFact_redist*elementDiameter.data()[eN]/degree_polynomial;
1824  // compute (smoothed) heaviside functions //
1825  double Hu = heaviside(u);
1826  double Huh = heaviside(uh);
1827  // compute cell metrics //
1828  cell_I_err += fabs(Hu - Huh)*dV;
1829  cell_V += Huh*dV;
1830  cell_V0 += Hu*dV;
1831 
1832  double norm2_grad_uh = 0.;
1833  for (int I=0; I<nSpace; I++)
1834  norm2_grad_uh += grad_uh[I]*grad_uh[I];
1835  cell_D_err += std::pow(std::sqrt(norm2_grad_uh) - 1, 2.)*dV;
1836  }
1837  global_V += cell_V;
1838  global_V0 += cell_V0;
1839  // metrics //
1840  global_I_err += cell_I_err;
1841  global_D_err += cell_D_err;
1842  }//elements
1843  global_V_err = fabs(global_V0 - global_V)/global_V0;
1844  global_D_err *= 0.5;
1845  return std::tuple<double, double, double>(global_I_err, global_V_err, global_D_err);
1846  }
1847 
1848  };//RDLS
1849  inline RDLS_base* newRDLS(int nSpaceIn,
1850  int nQuadraturePoints_elementIn,
1851  int nDOF_mesh_trial_elementIn,
1852  int nDOF_trial_elementIn,
1853  int nDOF_test_elementIn,
1854  int nQuadraturePoints_elementBoundaryIn,
1855  int CompKernelFlag)
1856  {
1857  if (nSpaceIn == 2)
1858  return proteus::chooseAndAllocateDiscretization2D<RDLS_base,RDLS,CompKernel>(nSpaceIn,
1859  nQuadraturePoints_elementIn,
1860  nDOF_mesh_trial_elementIn,
1861  nDOF_trial_elementIn,
1862  nDOF_test_elementIn,
1863  nQuadraturePoints_elementBoundaryIn,
1864  CompKernelFlag);
1865  else
1866  return proteus::chooseAndAllocateDiscretization<RDLS_base,RDLS,CompKernel>(nSpaceIn,
1867  nQuadraturePoints_elementIn,
1868  nDOF_mesh_trial_elementIn,
1869  nDOF_trial_elementIn,
1870  nDOF_test_elementIn,
1871  nQuadraturePoints_elementBoundaryIn,
1872  CompKernelFlag);
1873  }
1874 
1875 }//proteus
1876 
1877 #endif
proteus::RDLS_base
Definition: RDLS.h:28
proteus::RDLS::calculateJacobian_ellipticRedist
void calculateJacobian_ellipticRedist(arguments_dict &args, bool useExact)
Definition: RDLS.h:1416
proteus::RDLS::ck
CompKernelType ck
Definition: RDLS.h:51
SINGLE_POTENTIAL
#define SINGLE_POTENTIAL
Definition: RDLS.h:14
coeff
Double * coeff
Definition: Headers.h:42
proteus::RDLS::calculateMetricsAtEOS
std::tuple< double, double, double > calculateMetricsAtEOS(arguments_dict &args)
Definition: RDLS.h:1748
proteus::RDLS::evaluateCoefficients
void evaluateCoefficients(const double &eps, const double &u_levelSet, const double &u, const double grad_u[nSpace], double &m, double &dm, double &H, double dH[nSpace], double &r)
Definition: RDLS.h:59
proteus::RDLS::gfu
GeneralizedFunctions< nSpace, 2, nQuadraturePoints_element, nQuadraturePoints_elementBoundary > gfu
Definition: RDLS.h:52
proteus::RDLS::calculateResidual
void calculateResidual(arguments_dict &args, bool useExact)
Definition: RDLS.h:130
proteus::RDLS_base::calculateJacobian
virtual void calculateJacobian(arguments_dict &args, bool useExact)=0
proteus::RDLS_base::weighted_lumped_mass_matrix
std::valarray< double > weighted_lumped_mass_matrix
Definition: RDLS.h:30
ModelFactory.h
CompKernel.h
proteus::arguments_dict::scalar
T & scalar(const std::string &key)
proteus::RDLS_base::normalReconstruction
virtual void normalReconstruction(arguments_dict &args)=0
proteus::arguments_dict::array
xt::pyarray< T > & array(const std::string &key)
proteus::RDLS_base::calculateMetricsAtEOS
virtual std::tuple< double, double, double > calculateMetricsAtEOS(arguments_dict &args)=0
proteus::RDLS_base::~RDLS_base
virtual ~RDLS_base()
Definition: RDLS.h:31
proteus::RDLS_base::calculateResidual
virtual void calculateResidual(arguments_dict &args, bool useExact)=0
H
Double H
Definition: Headers.h:65
proteus::RDLS_base::calculateResidual_ellipticRedist
virtual void calculateResidual_ellipticRedist(arguments_dict &args, bool useExact)=0
proteus::RDLS::normalReconstruction
void normalReconstruction(arguments_dict &args)
Definition: RDLS.h:1624
equivalent_polynomials.h
proteus::RDLS::calculateJacobian
void calculateJacobian(arguments_dict &args, bool useExact)
Definition: RDLS.h:609
z
Double * z
Definition: Headers.h:49
proteus::RDLS_base::calculateJacobian_ellipticRedist
virtual void calculateJacobian_ellipticRedist(arguments_dict &args, bool useExact)=0
proteus::heaviside
double heaviside(const double &z)
Definition: CLSVOF.h:20
u
Double u
Definition: Headers.h:89
proteus::RDLS::calculateResidual_ellipticRedist
void calculateResidual_ellipticRedist(arguments_dict &args, bool useExact)
Definition: RDLS.h:1102
proteus::RDLS::nDOF_test_X_trial_element
const int nDOF_test_X_trial_element
Definition: RDLS.h:50
equivalent_polynomials::GeneralizedFunctions_mix
Definition: equivalent_polynomials.h:767
proteus::RDLS::gf
GeneralizedFunctions< nSpace, 2, nQuadraturePoints_element, nQuadraturePoints_elementBoundary > gf
Definition: RDLS.h:52
proteus
Definition: ADR.h:17
proteus::GeneralizedFunctions
equivalent_polynomials::GeneralizedFunctions_mix< nSpace, nP, nQ, nEBQ > GeneralizedFunctions
Definition: ADR.h:19
proteus::RDLS
Definition: RDLS.h:48
r
Double r
Definition: Headers.h:83
proteus::arguments_dict
Definition: ArgumentsDict.h:70
proteus::RDLS::calculateSubgridError_tau
void calculateSubgridError_tau(const double &elementDiameter, const double &dmt, const double dH[nSpace], double &cfl, double &tau)
Definition: RDLS.h:95
proteus::RDLS::calculateSubgridError_tau
void calculateSubgridError_tau(const double G[nSpace *nSpace], const double Ai[nSpace], double &tau_v, double &q_cfl)
Definition: RDLS.h:116
proteus::RDLS::RDLS
RDLS()
Definition: RDLS.h:53
proteus::newRDLS
RDLS_base * newRDLS(int nSpaceIn, int nQuadraturePoints_elementIn, int nDOF_mesh_trial_elementIn, int nDOF_trial_elementIn, int nDOF_test_elementIn, int nQuadraturePoints_elementBoundaryIn, int CompKernelFlag)
Definition: RDLS.h:1849
ArgumentsDict.h
Generated on Fri Jul 1 2022 11:17:40 for proteus by doxygen 1.8.20