proteus  1.8.1
C/C++/Fortran libraries
PresInc.h
Go to the documentation of this file.
1 #ifndef PRESINC_H
2 #define PRESINC_H
3 #include <cmath>
4 #include <iostream>
5 #include "CompKernel.h"
6 #include "ModelFactory.h"
7 #include "ArgumentsDict.h"
8 #include "xtensor-python/pyarray.hpp"
9 
10 namespace py = pybind11;
11 
12 namespace proteus
13 {
14  class cppPresInc_base
15  {
16  public:
17  virtual ~cppPresInc_base(){}
18  virtual void calculateResidual(arguments_dict& args)=0;
19  virtual void calculateJacobian(arguments_dict& args)=0;
20  };
21 
22  template<class CompKernelType,
23  int nSpace,
24  int nQuadraturePoints_element,
25  int nDOF_mesh_trial_element,
26  int nDOF_trial_element,
27  int nDOF_test_element,
28  int nQuadraturePoints_elementBoundary>
29  class cppPresInc : public cppPresInc_base
30  {
31  public:
32  const int nDOF_test_X_trial_element;
33  CompKernelType ck;
34  cppPresInc():
35  nDOF_test_X_trial_element(nDOF_test_element*nDOF_trial_element),
36  ck()
37  {}
38  inline
39  void evaluateCoefficients(const double& alphaBDF,
40  const double vf[nSpace],
41  const double vs[nSpace],
42  const double& vos,
43  const double& rhos_min,
44  const double& rhof_min,
45  double f[nSpace],
46  double& a)
47  {
48  for (int I=0;I<nSpace;I++)
49  f[I] = (1.0-vos)*vf[I] + vos*vs[I];
50  //a = 1.0/( vos*rhos_min*alphaBDF) + 1.0/( (1.0-vos)*rhof_min*alphaBDF );
51  a = 1.0/( vos*rhos_min*alphaBDF + (1.0-vos)*rhof_min*alphaBDF );
52  //a = (1.0-vos)/(rhof_min*alphaBDF) + vos*(rhos_min/alphaBDF);
53  //a = (1.0)/(rhof_min*alphaBDF) ;
54  //a = (1.0-vos)/(rhof_min*alphaBDF) + (vos)/(rhos_min*alphaBDF);
55  }
56 
57  inline
58  void exteriorNumericalAdvectiveFlux(const int& isFluxBoundary,
59  const double& bc_flux,
60  const double n[nSpace],
61  const double f[nSpace],
62  double& flux)
63  {
64  if (isFluxBoundary == 1)
65  flux = bc_flux;
66  else
67  {
68  flux = 0.0;
69  for (int I=0; I < nSpace; I++)
70  flux += n[I]*f[I];
71  }
72  }
73 
74  inline
75  void exteriorNumericalDiffusiveFlux(const int& isDOFBoundary,
76  const int& isFluxBoundary,
77  const double n[nSpace],
78  const double& a,
79  const double grad_potential[nSpace],
80  const double& u,
81  const double& bc_u,
82  const double& bc_flux,
83  const double& penalty,
84  double& flux)
85  {
86  if(isFluxBoundary == 1)
87  {
88  flux = bc_flux;
89  }
90  else if(isDOFBoundary == 1)
91  {
92  flux = 0.0;
93  for(int I=0;I<nSpace;I++)
94  flux-= a*grad_potential[I]*n[I];
95  flux += a*penalty*(u-bc_u);
96  }
97  else
98  {
99  std::cerr<<"PresInc.h: warning, diffusion term with no boundary condition set, setting diffusive flux to 0.0"<<std::endl;
100  flux = 0.0;
101  }
102  }
103 
104  inline
105  double ExteriorNumericalDiffusiveFluxJacobian(const int& isDOFBoundary,
106  const int& isFluxBoundary,
107  const double n[nSpace],
108  const double& a,
109  const double& v,
110  const double grad_v[nSpace],
111  const double& penalty)
112  {
113  double tmp=0.0;
114  if(isFluxBoundary==0 && isDOFBoundary==1)
115  {
116  for(int I=0;I<nSpace;I++)
117  tmp -= a*grad_v[I]*n[I];
118  tmp +=a*penalty*v;
119  }
120  return tmp;
121  }
122 
123  inline void calculateElementResidual(//element
124  double* mesh_trial_ref,
125  double* mesh_grad_trial_ref,
126  double* mesh_dof,
127  int* mesh_l2g,
128  double* dV_ref,
129  double* u_trial_ref,
130  double* u_grad_trial_ref,
131  double* u_test_ref,
132  double* u_grad_test_ref,
133  //element boundary
134  double* mesh_trial_trace_ref,
135  double* mesh_grad_trial_trace_ref,
136  double* dS_ref,
137  double* u_trial_trace_ref,
138  double* u_grad_trial_trace_ref,
139  double* u_test_trace_ref,
140  double* u_grad_test_trace_ref,
141  double* normal_ref,
142  double* boundaryJac_ref,
143  //physics
144  int nElements_global,
145  int* u_l2g,
146  double* u_dof,
147  double alphaBDF,
148  double* q_vf,
149  double* q_divU,
150  double* q_vs,
151  double* q_vos,
152  double rho_s,
153  double* q_rho_f,
154  double rho_s_min,
155  double rho_f_min,
156  double* ebqe_vf,
157  double* ebqe_vs,
158  double* ebqe_vos,
159  double* ebqe_rho_f,
160  double* q_u,
161  double* q_grad_u,
162  double* ebqe_u,
163  double* ebqe_grad_u,
164  int offset_u,
165  int stride_u,
166  double* elementResidual_u,
167  int nExteriorElementBoundaries_global,
168  int* exteriorElementBoundariesArray,
169  int* elementBoundaryElementsArray,
170  int* elementBoundaryLocalElementBoundariesArray,
171  double* element_u,
172  int eN,
173  double compatibility_condition,
174  int INTEGRATE_BY_PARTS_DIV_U,
175  double* q_a)
176  {
177  for (int i=0;i<nDOF_test_element;i++)
178  {
179  elementResidual_u[i]=0.0;
180  }//i
181  //loop over quadrature points and compute integrands
182  for (int k=0;k<nQuadraturePoints_element;k++)
183  {
184  //compute indeces and declare local storage
185  register int eN_k = eN*nQuadraturePoints_element+k,
186  eN_k_nSpace = eN_k*nSpace;
187  //eN_nDOF_trial_element = eN*nDOF_trial_element;
188  register double u=0.0,grad_u[nSpace],
189  a=0.0,
190  f[nSpace],
191  jac[nSpace*nSpace],
192  jacDet,
193  jacInv[nSpace*nSpace],
194  u_grad_trial[nDOF_trial_element*nSpace],
195  u_test_dV[nDOF_trial_element],
196  u_grad_test_dV[nDOF_test_element*nSpace],
197  dV,x,y,z,
198  G[nSpace*nSpace],G_dd_G,tr_G;
199  //
200  //compute solution and gradients at quadrature points
201  //
202  ck.calculateMapping_element(eN,
203  k,
204  mesh_dof,
205  mesh_l2g,
206  mesh_trial_ref,
207  mesh_grad_trial_ref,
208  jac,
209  jacDet,
210  jacInv,
211  x,y,z);
212  //get the physical integration weight
213  dV = fabs(jacDet)*dV_ref[k];
214  ck.calculateG(jacInv,G,G_dd_G,tr_G);
215  //get the trial function gradients
216  ck.gradTrialFromRef(&u_grad_trial_ref[k*nDOF_trial_element*nSpace],jacInv,u_grad_trial);
217  //get the solution
218  ck.valFromElementDOF(element_u,&u_trial_ref[k*nDOF_trial_element],u);
219  //get the solution gradients
220  ck.gradFromElementDOF(element_u,u_grad_trial,grad_u);
221  //precalculate test function products with integration weights
222  for (int j=0;j<nDOF_trial_element;j++)
223  {
224  u_test_dV[j] = u_test_ref[k*nDOF_trial_element+j]*dV;
225  for (int I=0;I<nSpace;I++)
226  {
227  u_grad_test_dV[j*nSpace+I] = u_grad_trial[j*nSpace+I]*dV;//cek warning won't work for Petrov-Galerkin
228  }
229  }
230  //
231  //calculate pde coefficients at quadrature points
232  //
233  evaluateCoefficients(alphaBDF,
234  &q_vf[eN_k_nSpace],
235  &q_vs[eN_k_nSpace],
236  q_vos[eN_k],
237  rho_s_min,
238  rho_f_min,
239  f,
240  a);
241  //
242  //update element residual
243  //
244  for(int i=0;i<nDOF_test_element;i++)
245  {
246  //register int eN_k_i=eN_k*nDOF_test_element+i;
247  //register int eN_k_i_nSpace = eN_k_i*nSpace;
248  register int i_nSpace=i*nSpace;
249  elementResidual_u[i] +=
250  (INTEGRATE_BY_PARTS_DIV_U == 1 ? ck.Advection_weak(f,&u_grad_test_dV[i_nSpace]) : q_divU[eN_k]*u_test_dV[i])
251  //+ compatibility_condition*u_test_dV[i] // mql: to make the system solvable if int(div(u))!=0
252  + ck.NumericalDiffusion(a,grad_u,&u_grad_test_dV[i_nSpace]);
253  }//i
254  //
255  //save momentum for time history and velocity for subgrid error
256  //save solution for other models
257  //
258  q_a[eN_k] = a;
259  q_u[eN_k] = u;
260 
261  for (int I=0;I<nSpace;I++)
262  q_grad_u[eN_k_nSpace+I] = grad_u[I];
263  }
264  }
265 
266  void calculateResidual(arguments_dict& args)
267  {
268  xt::pyarray<double>& mesh_trial_ref = args.array<double>("mesh_trial_ref");
269  xt::pyarray<double>& mesh_grad_trial_ref = args.array<double>("mesh_grad_trial_ref");
270  xt::pyarray<double>& mesh_dof = args.array<double>("mesh_dof");
271  xt::pyarray<int>& mesh_l2g = args.array<int>("mesh_l2g");
272  xt::pyarray<double>& dV_ref = args.array<double>("dV_ref");
273  xt::pyarray<double>& u_trial_ref = args.array<double>("u_trial_ref");
274  xt::pyarray<double>& u_grad_trial_ref = args.array<double>("u_grad_trial_ref");
275  xt::pyarray<double>& u_test_ref = args.array<double>("u_test_ref");
276  xt::pyarray<double>& u_grad_test_ref = args.array<double>("u_grad_test_ref");
277  xt::pyarray<double>& mesh_trial_trace_ref = args.array<double>("mesh_trial_trace_ref");
278  xt::pyarray<double>& mesh_grad_trial_trace_ref = args.array<double>("mesh_grad_trial_trace_ref");
279  xt::pyarray<double>& dS_ref = args.array<double>("dS_ref");
280  xt::pyarray<double>& u_trial_trace_ref = args.array<double>("u_trial_trace_ref");
281  xt::pyarray<double>& u_grad_trial_trace_ref = args.array<double>("u_grad_trial_trace_ref");
282  xt::pyarray<double>& u_test_trace_ref = args.array<double>("u_test_trace_ref");
283  xt::pyarray<double>& u_grad_test_trace_ref = args.array<double>("u_grad_test_trace_ref");
284  xt::pyarray<double>& normal_ref = args.array<double>("normal_ref");
285  xt::pyarray<double>& boundaryJac_ref = args.array<double>("boundaryJac_ref");
286  int nElements_global = args.scalar<int>("nElements_global");
287  xt::pyarray<int>& isDOFBoundary = args.array<int>("isDOFBoundary");
288  xt::pyarray<int>& isFluxBoundary = args.array<int>("isFluxBoundary");
289  xt::pyarray<int>& u_l2g = args.array<int>("u_l2g");
290  xt::pyarray<double>& u_dof = args.array<double>("u_dof");
291  double alphaBDF = args.scalar<double>("alphaBDF");
292  xt::pyarray<double>& q_vf = args.array<double>("q_vf");
293  xt::pyarray<double>& q_divU = args.array<double>("q_divU");
294  xt::pyarray<double>& q_vs = args.array<double>("q_vs");
295  xt::pyarray<double>& q_vos = args.array<double>("q_vos");
296  double rho_s = args.scalar<double>("rho_s");
297  xt::pyarray<double>& q_rho_f = args.array<double>("q_rho_f");
298  double rho_s_min = args.scalar<double>("rho_s_min");
299  double rho_f_min = args.scalar<double>("rho_f_min");
300  xt::pyarray<double>& ebqe_vf = args.array<double>("ebqe_vf");
301  xt::pyarray<double>& ebqe_vs = args.array<double>("ebqe_vs");
302  xt::pyarray<double>& ebqe_vos = args.array<double>("ebqe_vos");
303  xt::pyarray<double>& ebqe_rho_f = args.array<double>("ebqe_rho_f");
304  xt::pyarray<double>& q_u = args.array<double>("q_u");
305  xt::pyarray<double>& q_grad_u = args.array<double>("q_grad_u");
306  xt::pyarray<double>& ebqe_u = args.array<double>("ebqe_u");
307  xt::pyarray<double>& ebqe_grad_u = args.array<double>("ebqe_grad_u");
308  xt::pyarray<double>& ebqe_bc_u_ext = args.array<double>("ebqe_bc_u_ext");
309  xt::pyarray<double>& ebqe_adv_flux = args.array<double>("ebqe_adv_flux");
310  xt::pyarray<double>& ebqe_diff_flux = args.array<double>("ebqe_diff_flux");
311  xt::pyarray<double>& bc_adv_flux = args.array<double>("bc_adv_flux");
312  xt::pyarray<double>& bc_diff_flux = args.array<double>("bc_diff_flux");
313  int offset_u = args.scalar<int>("offset_u");
314  int stride_u = args.scalar<int>("stride_u");
315  xt::pyarray<double>& globalResidual = args.array<double>("globalResidual");
316  int nExteriorElementBoundaries_global = args.scalar<int>("nExteriorElementBoundaries_global");
317  xt::pyarray<int>& exteriorElementBoundariesArray = args.array<int>("exteriorElementBoundariesArray");
318  xt::pyarray<int>& elementBoundaryElementsArray = args.array<int>("elementBoundaryElementsArray");
319  xt::pyarray<int>& elementBoundaryLocalElementBoundariesArray = args.array<int>("elementBoundaryLocalElementBoundariesArray");
320  int INTEGRATE_BY_PARTS_DIV_U = args.scalar<int>("INTEGRATE_BY_PARTS_DIV_U");
321  xt::pyarray<double>& q_a = args.array<double>("q_a");
322  xt::pyarray<double>& ebqe_a = args.array<double>("ebqe_a");
323  double compatibility_condition=0.;
324  /* // COMPUTE COMPATIBILITY CONSTANT */
325  /* // mql: Modify the rhs (by adding a constant) so that the Poission system is solvable (assume diffusive flux = 0). */
326  /* // Note that this is equivalent to consider a (not known) diffusive flux != 0 s.t. the system is solvable. */
327  /* for(int eN=0;eN<nElements_global;eN++) */
328  /* { */
329  /* for (int k=0;k<nQuadraturePoints_element;k++) */
330  /* { */
331  /* register int eN_k = eN*nQuadraturePoints_element+k; */
332  /* register double */
333  /* jac[nSpace*nSpace], */
334  /* jacDet, */
335  /* jacInv[nSpace*nSpace], */
336  /* dV,x,y,z; */
337  /* ck.calculateMapping_element(eN, */
338  /* k, */
339  /* mesh_dof.data(), */
340  /* mesh_l2g.data(), */
341  /* mesh_trial_ref.data(), */
342  /* mesh_grad_trial_ref.data(), */
343  /* jac, */
344  /* jacDet, */
345  /* jacInv, */
346  /* x,y,z); */
347  /* //get the physical integration weight */
348  /* dV = fabs(jacDet)*dV_ref.data()[k]; */
349  /* compatibility_condition -= q_divU.data()[eN_k]*dV; */
350  /* } */
351  /* } */
352  //
353  //loop over elements to compute volume integrals and load them into element and global residual
354  //
355  //eN is the element index
356  //eN_k is the quadrature point index for a scalar
357  //eN_k_nSpace is the quadrature point index for a vector
358  //eN_i is the element test function index
359  //eN_j is the element trial function index
360  //eN_k_j is the quadrature point index for a trial function
361  //eN_k_i is the quadrature point index for a trial function
362  for(int eN=0;eN<nElements_global;eN++)
363  {
364  //declare local storage for element residual and initialize
365  register double elementResidual_u[nDOF_test_element],element_u[nDOF_trial_element];
366  for (int i=0;i<nDOF_test_element;i++)
367  {
368  register int eN_i=eN*nDOF_test_element+i;
369  element_u[i] = u_dof.data()[u_l2g.data()[eN_i]];
370  }//i
371  calculateElementResidual(mesh_trial_ref.data(),
372  mesh_grad_trial_ref.data(),
373  mesh_dof.data(),
374  mesh_l2g.data(),
375  dV_ref.data(),
376  u_trial_ref.data(),
377  u_grad_trial_ref.data(),
378  u_test_ref.data(),
379  u_grad_test_ref.data(),
380  mesh_trial_trace_ref.data(),
381  mesh_grad_trial_trace_ref.data(),
382  dS_ref.data(),
383  u_trial_trace_ref.data(),
384  u_grad_trial_trace_ref.data(),
385  u_test_trace_ref.data(),
386  u_grad_test_trace_ref.data(),
387  normal_ref.data(),
388  boundaryJac_ref.data(),
389  nElements_global,
390  u_l2g.data(),
391  u_dof.data(),
392  alphaBDF,
393  q_vf.data(),
394  q_divU.data(),
395  q_vs.data(),
396  q_vos.data(),
397  rho_s,
398  q_rho_f.data(),
399  rho_s_min,
400  rho_f_min,
401  ebqe_vf.data(),
402  ebqe_vs.data(),
403  ebqe_vos.data(),
404  ebqe_rho_f.data(),
405  q_u.data(),
406  q_grad_u.data(),
407  ebqe_u.data(),
408  ebqe_grad_u.data(),
409  offset_u,
410  stride_u,
411  elementResidual_u,
412  nExteriorElementBoundaries_global,
413  exteriorElementBoundariesArray.data(),
414  elementBoundaryElementsArray.data(),
415  elementBoundaryLocalElementBoundariesArray.data(),
416  element_u,
417  eN,
418  compatibility_condition,
419  INTEGRATE_BY_PARTS_DIV_U,
420  q_a.data());
421  //
422  //load element into global residual and save element residual
423  //
424  for(int i=0;i<nDOF_test_element;i++)
425  {
426  register int eN_i=eN*nDOF_test_element+i;
427  globalResidual.data()[offset_u+stride_u*u_l2g.data()[eN_i]]+=elementResidual_u[i];
428  }//i
429  }//elements
430  //
431  //loop over exterior element boundaries to calculate levelset gradient
432  //
433  //ebNE is the Exterior element boundary INdex
434  //ebN is the element boundary INdex
435  //eN is the element index
436  for (int ebNE = 0; ebNE < nExteriorElementBoundaries_global; ebNE++)
437  {
438  register int ebN = exteriorElementBoundariesArray.data()[ebNE],
439  eN = elementBoundaryElementsArray.data()[ebN*2+0],
440  ebN_local = elementBoundaryLocalElementBoundariesArray.data()[ebN*2+0];
441  //eN_nDOF_trial_element = eN*nDOF_trial_element;
442  register double elementResidual_u[nDOF_test_element];
443  double element_u[nDOF_trial_element];
444  for (int i=0;i<nDOF_test_element;i++)
445  {
446  register int eN_i=eN*nDOF_test_element+i;
447  element_u[i] = u_dof.data()[u_l2g.data()[eN_i]];
448  elementResidual_u[i] = 0.0;
449  }//i
450  for (int kb=0;kb<nQuadraturePoints_elementBoundary;kb++)
451  {
452  register int ebNE_kb = ebNE*nQuadraturePoints_elementBoundary+kb,
453  ebNE_kb_nSpace = ebNE_kb*nSpace,
454  ebN_local_kb = ebN_local*nQuadraturePoints_elementBoundary+kb,
455  ebN_local_kb_nSpace = ebN_local_kb*nSpace;
456  register double h_b=0.0,
457  penalty=0.0,
458  u_ext=0.0,
459  bc_u_ext=0.0,
460  adv_flux_ext=0.0,
461  diff_flux_ext=0.0,
462  a_ext=0.0,
463  f_ext[nSpace],
464  grad_u_ext[nSpace],
465  jac_ext[nSpace*nSpace],
466  jacDet_ext,
467  jacInv_ext[nSpace*nSpace],
468  boundaryJac[nSpace*(nSpace-1)],
469  metricTensor[(nSpace-1)*(nSpace-1)],
470  metricTensorDetSqrt,
471  dS,
472  u_test_dS[nDOF_test_element],
473  u_grad_trial_trace[nDOF_trial_element*nSpace],
474  u_grad_test_dS[nDOF_test_element*nSpace],
475  normal[nSpace],x_ext,y_ext,z_ext,
476  G[nSpace*nSpace],G_dd_G,tr_G;
477  //
478  //calculate the solution and gradients at quadrature points
479  //
480  ck.calculateMapping_elementBoundary(eN,
481  ebN_local,
482  kb,
483  ebN_local_kb,
484  mesh_dof.data(),
485  mesh_l2g.data(),
486  mesh_trial_trace_ref.data(),
487  mesh_grad_trial_trace_ref.data(),
488  boundaryJac_ref.data(),
489  jac_ext,
490  jacDet_ext,
491  jacInv_ext,
492  boundaryJac,
493  metricTensor,
494  metricTensorDetSqrt,
495  normal_ref.data(),
496  normal,
497  x_ext,y_ext,z_ext);
498  dS = metricTensorDetSqrt*dS_ref.data()[kb];
499  //get the metric tensor
500  //cek todo use symmetry
501  ck.calculateG(jacInv_ext,G,G_dd_G,tr_G);
502  ck.calculateGScale(G,normal,h_b);
503  penalty = 10.0/h_b;
504  //compute shape and solution information
505  //shape
506  ck.gradTrialFromRef(&u_grad_trial_trace_ref.data()[ebN_local_kb_nSpace*nDOF_trial_element],jacInv_ext,u_grad_trial_trace);
507  //solution and gradients
508  ck.valFromElementDOF(element_u,&u_trial_trace_ref.data()[ebN_local_kb*nDOF_test_element],u_ext);
509  ck.gradFromElementDOF(element_u,u_grad_trial_trace,grad_u_ext);
510  //precalculate test function products with integration weights
511  for (int j=0;j<nDOF_trial_element;j++)
512  {
513  u_test_dS[j] = u_test_trace_ref.data()[ebN_local_kb*nDOF_test_element+j]*dS;
514  for (int I=0;I<nSpace;I++)
515  u_grad_test_dS[j*nSpace+I] = u_grad_trial_trace[j*nSpace+I]*dS;//cek hack, using trial
516  }
517  //
518  //load the boundary values
519  //
520  bc_u_ext = isDOFBoundary.data()[ebNE_kb]*ebqe_bc_u_ext.data()[ebNE_kb]+(1-isDOFBoundary.data()[ebNE_kb])*u_ext;
521  //
522  //calculate the pde coefficients using the solution and the boundary values for the solution
523  //
524 
525  evaluateCoefficients(alphaBDF,
526  &ebqe_vf.data()[ebNE_kb_nSpace],
527  &ebqe_vs.data()[ebNE_kb_nSpace],
528  q_vos.data()[ebNE_kb],
529  rho_s_min,
530  rho_f_min,
531  f_ext,
532  a_ext);
533  ebqe_u.data()[ebNE_kb] = u_ext;
534  for (int I=0;I<nSpace;I++)
535  ebqe_grad_u.data()[ebNE_kb_nSpace+I] = grad_u_ext[I];
536  //
537  //calculate the numerical fluxes
538  //
539  exteriorNumericalAdvectiveFlux(isFluxBoundary.data()[ebNE_kb],
540  bc_adv_flux.data()[ebNE_kb],
541  normal,
542  f_ext,
543  adv_flux_ext); //=f.normal = [(1-vos)*vf + vos*vs].normal
544  exteriorNumericalDiffusiveFlux(isDOFBoundary.data()[ebNE_kb],
545  isFluxBoundary.data()[ebNE_kb],
546  normal,
547  a_ext,
548  grad_u_ext,
549  u_ext,
550  bc_u_ext,
551  bc_diff_flux.data()[ebNE_kb],
552  penalty,
553  diff_flux_ext);
554  if (isDOFBoundary.data()[ebNE_kb] != 1)
555  diff_flux_ext = 0.0; // mql: don't consider diffusive flux unless Dirichlet BC
556  //if(isFluxBoundary.data()[ebNE_kb] == 1)
557  //{
558  // adv_flux_ext = 0.0;
559  // diff_flux_ext = 0.0;
560  //}
561  ebqe_a.data()[ebNE_kb] = a_ext;
562  ebqe_adv_flux.data()[ebNE_kb] = adv_flux_ext;
563  ebqe_diff_flux.data()[ebNE_kb] = diff_flux_ext;
564  //
565  //update residuals
566  //
567  for (int i=0;i<nDOF_test_element;i++)
568  {
569  elementResidual_u[i] +=
570  (INTEGRATE_BY_PARTS_DIV_U == 1 ? ck.ExteriorElementBoundaryFlux(adv_flux_ext,u_test_dS[i]) : 0.)
571  + ck.ExteriorElementBoundaryFlux(diff_flux_ext,u_test_dS[i]) // mql: just != 0 if Dirichlet BC
572  + ck.ExteriorElementBoundaryScalarDiffusionAdjoint(isDOFBoundary.data()[ebNE_kb],
573  0, // mql: if Dirichlet BCs, then the flux BCs don't matter
574  //isFluxBoundary.data()[ebNE_kb],
575  1.0,
576  u_ext,
577  bc_u_ext,
578  normal,
579  a_ext,
580  &u_grad_test_dS[i*nSpace]);
581  }//i
582  }//kb
583  //
584  //update the element and global residual storage
585  //
586  for (int i=0;i<nDOF_test_element;i++)
587  {
588  int eN_i = eN*nDOF_test_element+i;
589  globalResidual.data()[offset_u+stride_u*u_l2g.data()[eN_i]] += elementResidual_u[i];
590  }//i
591  }//ebNE
592 
593  }
594 
595  inline void calculateElementJacobian(//element
596  double* mesh_trial_ref,
597  double* mesh_grad_trial_ref,
598  double* mesh_dof,
599  int* mesh_l2g,
600  double* dV_ref,
601  double* u_trial_ref,
602  double* u_grad_trial_ref,
603  double* u_test_ref,
604  double* u_grad_test_ref,
605  //element boundary
606  double* mesh_trial_trace_ref,
607  double* mesh_grad_trial_trace_ref,
608  double* dS_ref,
609  double* u_trial_trace_ref,
610  double* u_grad_trial_trace_ref,
611  double* u_test_trace_ref,
612  double* u_grad_test_trace_ref,
613  double* normal_ref,
614  double* boundaryJac_ref,
615  //physics
616  int nElements_global,
617  int* u_l2g,
618  double* u_dof,
619  double alphaBDF,
620  double* q_vf,
621  double* q_vs,
622  double* q_vos,
623  double rho_s,
624  double* q_rho_f,
625  double rho_s_min,
626  double rho_f_min,
627  double* elementJacobian_u_u,
628  double* element_u,
629  int eN)
630  {
631  for (int i=0;i<nDOF_test_element;i++)
632  for (int j=0;j<nDOF_trial_element;j++)
633  {
634  elementJacobian_u_u[i*nDOF_trial_element+j]=0.0;
635  }
636  for (int k=0;k<nQuadraturePoints_element;k++)
637  {
638  int eN_k = eN*nQuadraturePoints_element+k, //index to a scalar at a quadrature point
639  eN_k_nSpace = eN_k*nSpace;
640  //eN_nDOF_trial_element = eN*nDOF_trial_element; //index to a vector at a quadrature point
641 
642  //declare local storage
643  register double u=0.0,
644  grad_u[nSpace],
645  f[nSpace],
646  a=0.0,
647  jac[nSpace*nSpace],
648  jacDet,
649  jacInv[nSpace*nSpace],
650  u_grad_trial[nDOF_trial_element*nSpace],
651  dV,
652  u_grad_test_dV[nDOF_test_element*nSpace],
653  x,y,z,
654  G[nSpace*nSpace],G_dd_G,tr_G;
655  //
656  //calculate solution and gradients at quadrature points
657  //
658  ck.calculateMapping_element(eN,
659  k,
660  mesh_dof,
661  mesh_l2g,
662  mesh_trial_ref,
663  mesh_grad_trial_ref,
664  jac,
665  jacDet,
666  jacInv,
667  x,y,z);
668  //get the physical integration weight
669  dV = fabs(jacDet)*dV_ref[k];
670  ck.calculateG(jacInv,G,G_dd_G,tr_G);
671  //get the trial function gradients
672  ck.gradTrialFromRef(&u_grad_trial_ref[k*nDOF_trial_element*nSpace],jacInv,u_grad_trial);
673  //get the solution
674  ck.valFromElementDOF(element_u,&u_trial_ref[k*nDOF_trial_element],u);
675  //get the solution gradients
676  ck.gradFromElementDOF(element_u,u_grad_trial,grad_u);
677  //precalculate test function products with integration weights
678  for (int j=0;j<nDOF_trial_element;j++)
679  {
680  for (int I=0;I<nSpace;I++)
681  {
682  u_grad_test_dV[j*nSpace+I] = u_grad_trial[j*nSpace+I]*dV;//cek warning won't work for Petrov-Galerkin
683  }
684  }
685  //
686  //calculate pde coefficients and derivatives at quadrature points
687  //
688  evaluateCoefficients(alphaBDF,
689  &q_vf[eN_k_nSpace],
690  &q_vs[eN_k_nSpace],
691  q_vos[eN_k],
692  rho_s_min,
693  rho_f_min,
694  f,
695  a);
696  for(int i=0;i<nDOF_test_element;i++)
697  {
698  //int eN_k_i=eN_k*nDOF_test_element+i;
699  //int eN_k_i_nSpace=eN_k_i*nSpace;
700  int i_nSpace=i*nSpace;
701  for(int j=0;j<nDOF_trial_element;j++)
702  {
703  //int eN_k_j=eN_k*nDOF_trial_element+j;
704  //int eN_k_j_nSpace = eN_k_j*nSpace;
705  int j_nSpace = j*nSpace;
706  elementJacobian_u_u[i*nDOF_trial_element+j] +=
707  ck.NumericalDiffusionJacobian(a,&u_grad_trial[j_nSpace],&u_grad_test_dV[i_nSpace]);
708  }//j
709  }//i
710  }//k
711  }
712 
713  void calculateJacobian(arguments_dict& args)
714  {
715  xt::pyarray<double>& mesh_trial_ref = args.array<double>("mesh_trial_ref");
716  xt::pyarray<double>& mesh_grad_trial_ref = args.array<double>("mesh_grad_trial_ref");
717  xt::pyarray<double>& mesh_dof = args.array<double>("mesh_dof");
718  xt::pyarray<int>& mesh_l2g = args.array<int>("mesh_l2g");
719  xt::pyarray<double>& dV_ref = args.array<double>("dV_ref");
720  xt::pyarray<double>& u_trial_ref = args.array<double>("u_trial_ref");
721  xt::pyarray<double>& u_grad_trial_ref = args.array<double>("u_grad_trial_ref");
722  xt::pyarray<double>& u_test_ref = args.array<double>("u_test_ref");
723  xt::pyarray<double>& u_grad_test_ref = args.array<double>("u_grad_test_ref");
724  xt::pyarray<double>& mesh_trial_trace_ref = args.array<double>("mesh_trial_trace_ref");
725  xt::pyarray<double>& mesh_grad_trial_trace_ref = args.array<double>("mesh_grad_trial_trace_ref");
726  xt::pyarray<double>& dS_ref = args.array<double>("dS_ref");
727  xt::pyarray<double>& u_trial_trace_ref = args.array<double>("u_trial_trace_ref");
728  xt::pyarray<double>& u_grad_trial_trace_ref = args.array<double>("u_grad_trial_trace_ref");
729  xt::pyarray<double>& u_test_trace_ref = args.array<double>("u_test_trace_ref");
730  xt::pyarray<double>& u_grad_test_trace_ref = args.array<double>("u_grad_test_trace_ref");
731  xt::pyarray<double>& normal_ref = args.array<double>("normal_ref");
732  xt::pyarray<double>& boundaryJac_ref = args.array<double>("boundaryJac_ref");
733  int nElements_global = args.scalar<int>("nElements_global");
734  xt::pyarray<int>& isDOFBoundary = args.array<int>("isDOFBoundary");
735  xt::pyarray<int>& isFluxBoundary = args.array<int>("isFluxBoundary");
736  xt::pyarray<int>& u_l2g = args.array<int>("u_l2g");
737  xt::pyarray<double>& u_dof = args.array<double>("u_dof");
738  double alphaBDF = args.scalar<double>("alphaBDF");
739  xt::pyarray<double>& q_vf = args.array<double>("q_vf");
740  xt::pyarray<double>& q_vs = args.array<double>("q_vs");
741  xt::pyarray<double>& q_vos = args.array<double>("q_vos");
742  double rho_s = args.scalar<double>("rho_s");
743  xt::pyarray<double>& q_rho_f = args.array<double>("q_rho_f");
744  double rho_s_min = args.scalar<double>("rho_s_min");
745  double rho_f_min = args.scalar<double>("rho_f_min");
746  xt::pyarray<double>& ebqe_vf = args.array<double>("ebqe_vf");
747  xt::pyarray<double>& ebqe_vs = args.array<double>("ebqe_vs");
748  xt::pyarray<double>& ebqe_vos = args.array<double>("ebqe_vos");
749  xt::pyarray<double>& ebqe_rho_f = args.array<double>("ebqe_rho_f");
750  xt::pyarray<int>& csrRowIndeces_u_u = args.array<int>("csrRowIndeces_u_u");
751  xt::pyarray<int>& csrColumnOffsets_u_u = args.array<int>("csrColumnOffsets_u_u");
752  xt::pyarray<double>& globalJacobian = args.array<double>("globalJacobian");
753  int nExteriorElementBoundaries_global = args.scalar<int>("nExteriorElementBoundaries_global");
754  xt::pyarray<int>& exteriorElementBoundariesArray = args.array<int>("exteriorElementBoundariesArray");
755  xt::pyarray<int>& elementBoundaryElementsArray = args.array<int>("elementBoundaryElementsArray");
756  xt::pyarray<int>& elementBoundaryLocalElementBoundariesArray = args.array<int>("elementBoundaryLocalElementBoundariesArray");
757  xt::pyarray<int>& csrColumnOffsets_eb_u_u = args.array<int>("csrColumnOffsets_eb_u_u");
758  //
759  //loop over elements to compute volume integrals and load them into the element Jacobians and global Jacobian
760  //
761  for(int eN=0;eN<nElements_global;eN++)
762  {
763  register double elementJacobian_u_u[nDOF_test_element*nDOF_trial_element],element_u[nDOF_trial_element];
764  for (int j=0;j<nDOF_trial_element;j++)
765  {
766  register int eN_j = eN*nDOF_trial_element+j;
767  element_u[j] = u_dof.data()[u_l2g.data()[eN_j]];
768  }
769  calculateElementJacobian(mesh_trial_ref.data(),
770  mesh_grad_trial_ref.data(),
771  mesh_dof.data(),
772  mesh_l2g.data(),
773  dV_ref.data(),
774  u_trial_ref.data(),
775  u_grad_trial_ref.data(),
776  u_test_ref.data(),
777  u_grad_test_ref.data(),
778  mesh_trial_trace_ref.data(),
779  mesh_grad_trial_trace_ref.data(),
780  dS_ref.data(),
781  u_trial_trace_ref.data(),
782  u_grad_trial_trace_ref.data(),
783  u_test_trace_ref.data(),
784  u_grad_test_trace_ref.data(),
785  normal_ref.data(),
786  boundaryJac_ref.data(),
787  nElements_global,
788  u_l2g.data(),
789  u_dof.data(),
790  alphaBDF,
791  q_vf.data(),
792  q_vs.data(),
793  q_vos.data(),
794  rho_s,
795  q_rho_f.data(),
796  rho_s_min,
797  rho_f_min,
798  elementJacobian_u_u,
799  element_u,
800  eN);
801  //
802  //load into element Jacobian into global Jacobian
803  //
804  for (int i=0;i<nDOF_test_element;i++)
805  {
806  int eN_i = eN*nDOF_test_element+i;
807  for (int j=0;j<nDOF_trial_element;j++)
808  {
809  int eN_i_j = eN_i*nDOF_trial_element+j;
810  globalJacobian.data()[csrRowIndeces_u_u.data()[eN_i] + csrColumnOffsets_u_u.data()[eN_i_j]] += elementJacobian_u_u[i*nDOF_trial_element+j];
811  }//j
812  }//i
813  }//elements
814  //
815  //loop over exterior element boundaries to compute the surface integrals and load them into the global Jacobian
816  //
817  for (int ebNE = 0; ebNE < nExteriorElementBoundaries_global; ebNE++)
818  {
819  register int ebN = exteriorElementBoundariesArray.data()[ebNE];
820  register int eN = elementBoundaryElementsArray.data()[ebN*2+0],
821  ebN_local = elementBoundaryLocalElementBoundariesArray.data()[ebN*2+0],
822  eN_nDOF_trial_element = eN*nDOF_trial_element;
823  for (int kb=0;kb<nQuadraturePoints_elementBoundary;kb++)
824  {
825  register int ebNE_kb = ebNE*nQuadraturePoints_elementBoundary+kb,
826  ebNE_kb_nSpace = ebNE_kb*nSpace,
827  ebN_local_kb = ebN_local*nQuadraturePoints_elementBoundary+kb,
828  ebN_local_kb_nSpace = ebN_local_kb*nSpace;
829  register double h_b=0.0,
830  u_ext=0.0,
831  grad_u_ext[nSpace],
832  a_ext=0.0,
833  f_ext[nSpace],
834  jac_ext[nSpace*nSpace],
835  jacDet_ext,
836  jacInv_ext[nSpace*nSpace],
837  boundaryJac[nSpace*(nSpace-1)],
838  metricTensor[(nSpace-1)*(nSpace-1)],
839  metricTensorDetSqrt,
840  dS,
841  u_test_dS[nDOF_test_element],
842  u_grad_trial_trace[nDOF_trial_element*nSpace],
843  u_grad_test_dS[nDOF_test_element*nSpace],
844  normal[nSpace],x_ext,y_ext,z_ext,
845  penalty=0.0,
846  //
847  G[nSpace*nSpace],G_dd_G,tr_G;
848  //
849  //calculate the solution and gradients at quadrature points
850  //
851  u_ext=0.0;
852  for (int I=0;I<nSpace;I++)
853  {
854  grad_u_ext[I] = 0.0;
855  }
856  ck.calculateMapping_elementBoundary(eN,
857  ebN_local,
858  kb,
859  ebN_local_kb,
860  mesh_dof.data(),
861  mesh_l2g.data(),
862  mesh_trial_trace_ref.data(),
863  mesh_grad_trial_trace_ref.data(),
864  boundaryJac_ref.data(),
865  jac_ext,
866  jacDet_ext,
867  jacInv_ext,
868  boundaryJac,
869  metricTensor,
870  metricTensorDetSqrt,
871  normal_ref.data(),
872  normal,
873  x_ext,y_ext,z_ext);
874  dS = metricTensorDetSqrt*dS_ref.data()[kb];
875  ck.calculateG(jacInv_ext,G,G_dd_G,tr_G);
876  ck.calculateGScale(G,normal,h_b);
877  penalty=10.0/h_b;
878  //compute shape and solution information
879  //shape
880  ck.gradTrialFromRef(&u_grad_trial_trace_ref.data()[ebN_local_kb_nSpace*nDOF_trial_element],jacInv_ext,u_grad_trial_trace);
881  //solution and gradients
882  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);
883  ck.gradFromDOF(u_dof.data(),&u_l2g.data()[eN_nDOF_trial_element],u_grad_trial_trace,grad_u_ext);
884  //precalculate test function products with integration weights
885  for (int j=0;j<nDOF_trial_element;j++)
886  {
887  u_test_dS[j] = u_test_trace_ref.data()[ebN_local_kb*nDOF_test_element+j]*dS;
888  for (int I=0;I<nSpace;I++)
889  u_grad_test_dS[j*nSpace+I] = u_grad_trial_trace[j*nSpace+I]*dS;//cek hack, using trial
890  }
891  //
892  //calculate the internal and external trace of the pde coefficients
893  //
894  evaluateCoefficients(alphaBDF,
895  &ebqe_vf.data()[ebNE_kb_nSpace],
896  &ebqe_vs.data()[ebNE_kb_nSpace],
897  q_vos.data()[ebNE_kb],
898  rho_s_min,
899  rho_f_min,
900  f_ext,
901  a_ext);
902  //
903  //update the global Jacobian from the flux Jacobian
904  //
905  for (int i=0;i<nDOF_test_element;i++)
906  {
907  register int eN_i = eN*nDOF_test_element+i;
908  //register int ebNE_kb_i = ebNE_kb*nDOF_test_element+i;
909  for (int j=0;j<nDOF_trial_element;j++)
910  {
911  register int ebN_i_j = ebN*4*nDOF_test_X_trial_element + i*nDOF_trial_element + j,
912  ebN_local_kb_j=ebN_local_kb*nDOF_trial_element+j,
913  j_nSpace = j*nSpace;
914 
915  globalJacobian.data()[csrRowIndeces_u_u.data()[eN_i] + csrColumnOffsets_eb_u_u.data()[ebN_i_j]] +=
916  ExteriorNumericalDiffusiveFluxJacobian(isDOFBoundary.data()[ebNE_kb],
917  isFluxBoundary.data()[ebNE_kb],
918  normal,
919  a_ext,
920  u_trial_trace_ref.data()[ebN_local_kb_j],
921  &u_grad_trial_trace[j_nSpace],
922  penalty)*u_test_dS[i]
923  +
924  ck.ExteriorElementBoundaryScalarDiffusionAdjointJacobian
925  (isDOFBoundary.data()[ebNE_kb],
926  isFluxBoundary.data()[ebNE_kb],
927  1.0,
928  u_trial_trace_ref.data()[ebN_local_kb_j],
929  normal,
930  a_ext,
931  &u_grad_test_dS[i*nSpace]);
932  }//j
933  }//i
934  }//kb
935  }//ebNE
936  }//computeJacobian
937  };//cppPresInc
938 
939  inline cppPresInc_base* newPresInc(int nSpaceIn,
940  int nQuadraturePoints_elementIn,
941  int nDOF_mesh_trial_elementIn,
942  int nDOF_trial_elementIn,
943  int nDOF_test_elementIn,
944  int nQuadraturePoints_elementBoundaryIn,
945  int CompKernelFlag)
946  {
947  if (nSpaceIn == 2)
948  return proteus::chooseAndAllocateDiscretization2D<cppPresInc_base,cppPresInc,CompKernel>(nSpaceIn,
949  nQuadraturePoints_elementIn,
950  nDOF_mesh_trial_elementIn,
951  nDOF_trial_elementIn,
952  nDOF_test_elementIn,
953  nQuadraturePoints_elementBoundaryIn,
954  CompKernelFlag);
955  else
956  return proteus::chooseAndAllocateDiscretization<cppPresInc_base,cppPresInc,CompKernel>(nSpaceIn,
957  nQuadraturePoints_elementIn,
958  nDOF_mesh_trial_elementIn,
959  nDOF_trial_elementIn,
960  nDOF_test_elementIn,
961  nQuadraturePoints_elementBoundaryIn,
962  CompKernelFlag);
963  }
964 }//proteus
965 #endif
proteus::cppPresInc_base::~cppPresInc_base
virtual ~cppPresInc_base()
Definition: PresInc.h:17
proteus::newPresInc
cppPresInc_base * newPresInc(int nSpaceIn, int nQuadraturePoints_elementIn, int nDOF_mesh_trial_elementIn, int nDOF_trial_elementIn, int nDOF_test_elementIn, int nQuadraturePoints_elementBoundaryIn, int CompKernelFlag)
Definition: PresInc.h:939
proteus::cppPresInc::calculateElementResidual
void calculateElementResidual(double *mesh_trial_ref, double *mesh_grad_trial_ref, double *mesh_dof, int *mesh_l2g, double *dV_ref, double *u_trial_ref, double *u_grad_trial_ref, double *u_test_ref, double *u_grad_test_ref, double *mesh_trial_trace_ref, double *mesh_grad_trial_trace_ref, double *dS_ref, double *u_trial_trace_ref, double *u_grad_trial_trace_ref, double *u_test_trace_ref, double *u_grad_test_trace_ref, double *normal_ref, double *boundaryJac_ref, int nElements_global, int *u_l2g, double *u_dof, double alphaBDF, double *q_vf, double *q_divU, double *q_vs, double *q_vos, double rho_s, double *q_rho_f, double rho_s_min, double rho_f_min, double *ebqe_vf, double *ebqe_vs, double *ebqe_vos, double *ebqe_rho_f, double *q_u, double *q_grad_u, double *ebqe_u, double *ebqe_grad_u, int offset_u, int stride_u, double *elementResidual_u, int nExteriorElementBoundaries_global, int *exteriorElementBoundariesArray, int *elementBoundaryElementsArray, int *elementBoundaryLocalElementBoundariesArray, double *element_u, int eN, double compatibility_condition, int INTEGRATE_BY_PARTS_DIV_U, double *q_a)
Definition: PresInc.h:123
proteus::cppPresInc::calculateResidual
void calculateResidual(arguments_dict &args)
Definition: PresInc.h:266
proteus::cppPresInc_base
Definition: PresInc.h:15
proteus::cppPresInc::cppPresInc
cppPresInc()
Definition: PresInc.h:34
proteus::cppPresInc_base::calculateJacobian
virtual void calculateJacobian(arguments_dict &args)=0
proteus::cppPresInc::calculateJacobian
void calculateJacobian(arguments_dict &args)
Definition: PresInc.h:713
n
Int n
Definition: Headers.h:28
proteus::cppPresInc::calculateElementJacobian
void calculateElementJacobian(double *mesh_trial_ref, double *mesh_grad_trial_ref, double *mesh_dof, int *mesh_l2g, double *dV_ref, double *u_trial_ref, double *u_grad_trial_ref, double *u_test_ref, double *u_grad_test_ref, double *mesh_trial_trace_ref, double *mesh_grad_trial_trace_ref, double *dS_ref, double *u_trial_trace_ref, double *u_grad_trial_trace_ref, double *u_test_trace_ref, double *u_grad_test_trace_ref, double *normal_ref, double *boundaryJac_ref, int nElements_global, int *u_l2g, double *u_dof, double alphaBDF, double *q_vf, double *q_vs, double *q_vos, double rho_s, double *q_rho_f, double rho_s_min, double rho_f_min, double *elementJacobian_u_u, double *element_u, int eN)
Definition: PresInc.h:595
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)
proteus::cppPresInc::ExteriorNumericalDiffusiveFluxJacobian
double ExteriorNumericalDiffusiveFluxJacobian(const int &isDOFBoundary, const int &isFluxBoundary, const double n[nSpace], const double &a, const double &v, const double grad_v[nSpace], const double &penalty)
Definition: PresInc.h:105
proteus::cppPresInc
Definition: PresInc.h:30
v
Double v
Definition: Headers.h:95
proteus::cppPresInc::calculateResidual
void calculateResidual(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 > &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 > &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, xt::pyarray< int > &isDOFBoundary, xt::pyarray< int > &isFluxBoundary, xt::pyarray< int > &u_l2g, xt::pyarray< double > &u_dof, double alphaBDF, xt::pyarray< double > &q_vf, xt::pyarray< double > &q_divU, xt::pyarray< double > &q_vs, xt::pyarray< double > &q_vos, double rho_s, xt::pyarray< double > &q_rho_f, double rho_s_min, double rho_f_min, xt::pyarray< double > &ebqe_vf, xt::pyarray< double > &ebqe_vs, xt::pyarray< double > &ebqe_vos, xt::pyarray< double > &ebqe_rho_f, xt::pyarray< double > &q_u, xt::pyarray< double > &q_grad_u, xt::pyarray< double > &ebqe_u, xt::pyarray< double > &ebqe_grad_u, xt::pyarray< double > &ebqe_bc_u_ext, xt::pyarray< double > &ebqe_adv_flux, xt::pyarray< double > &ebqe_diff_flux, xt::pyarray< double > &bc_adv_flux, xt::pyarray< double > &bc_diff_flux, int offset_u, int stride_u, xt::pyarray< double > &globalResidual, int nExteriorElementBoundaries_global, xt::pyarray< int > &exteriorElementBoundariesArray, xt::pyarray< int > &elementBoundaryElementsArray, xt::pyarray< int > &elementBoundaryLocalElementBoundariesArray, int INTEGRATE_BY_PARTS_DIV_U, xt::pyarray< double > &q_a, xt::pyarray< double > &ebqe_a)
Definition: PresInc.h:367
z
Double * z
Definition: Headers.h:49
u
Double u
Definition: Headers.h:89
proteus::cppPresInc::exteriorNumericalAdvectiveFlux
void exteriorNumericalAdvectiveFlux(const int &isFluxBoundary, const double &bc_flux, const double n[nSpace], const double f[nSpace], double &flux)
Definition: PresInc.h:58
proteus::cppPresInc::evaluateCoefficients
void evaluateCoefficients(const double &alphaBDF, const double vf[nSpace], const double vs[nSpace], const double &vos, const double &rhos_min, const double &rhof_min, double f[nSpace], double &a)
Definition: PresInc.h:39
proteus::cppPresInc::ck
CompKernelType ck
Definition: PresInc.h:33
proteus
Definition: ADR.h:17
proteus::f
double f(const double &g, const double &h, const double &hZ)
Definition: SW2DCV.h:58
proteus::arguments_dict
Definition: ArgumentsDict.h:70
proteus::cppPresInc_base::calculateResidual
virtual void calculateResidual(arguments_dict &args)=0
proteus::cppPresInc::nDOF_test_X_trial_element
const int nDOF_test_X_trial_element
Definition: PresInc.h:32
proteus::cppPresInc::exteriorNumericalDiffusiveFlux
void exteriorNumericalDiffusiveFlux(const int &isDOFBoundary, const int &isFluxBoundary, const double n[nSpace], const double &a, const double grad_potential[nSpace], const double &u, const double &bc_u, const double &bc_flux, const double &penalty, double &flux)
Definition: PresInc.h:75
ArgumentsDict.h
Generated on Fri Jul 1 2022 11:17:39 for proteus by doxygen 1.8.20