proteus  1.8.1
C/C++/Fortran libraries
PresInit.h
Go to the documentation of this file.
1 #ifndef PRESINIT_H
2 #define PRESINIT_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 cppPresInit_base
15  {
16  public:
17  virtual ~cppPresInit_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 cppPresInit : public cppPresInit_base
30  {
31  public:
32  CompKernelType ck;
33  cppPresInit():ck()
34  {}
35  inline double smoothedHeaviside(double eps, double phi)
36  {
37  double H;
38  if (phi > eps)
39  H=1.0;
40  else if (phi < -eps)
41  H=0.0;
42  else if (phi==0.0)
43  H=0.5;
44  else
45  H = 0.5*(1.0 + phi/eps + sin(M_PI*phi/eps)/M_PI);
46  return H;
47  }
48 
49  inline double smoothedHeaviside_integral(double eps, double phi)
50  {
51  double HI;
52  if (phi > eps)
53  {
54  HI= phi - eps + 0.5*(eps + 0.5*eps*eps/eps - eps*cos(M_PI*eps/eps)/(M_PI*M_PI)) - 0.5*((-eps) + 0.5*(-eps)*(-eps)/eps - eps*cos(M_PI*(-eps)/eps)/(M_PI*M_PI));
55  }
56  else if (phi < -eps)
57  {
58  HI=0.0;
59  }
60  else
61  {
62  HI = 0.5*(phi + 0.5*phi*phi/eps - eps*cos(M_PI*phi/eps)/(M_PI*M_PI)) - 0.5*((-eps) + 0.5*(-eps)*(-eps)/eps - eps*cos(M_PI*(-eps)/eps)/(M_PI*M_PI));
63  }
64  return HI;
65  }
66 
67  inline double smoothedDirac(double eps, double phi)
68  {
69  double d;
70  if (phi > eps)
71  d=0.0;
72  else if (phi < -eps)
73  d=0.0;
74  else
75  d = 0.5*(1.0 + cos(M_PI*phi/eps))/eps;
76  return d;
77  }
78  inline
79  void evaluateCoefficients(const double& epsHeaviside,
80  const double& epsDirac,
81  const double& phi,
82  const double& H,
83  const double& u,
84  const double& porosity,
85  double& r,
86  double& dr)
87  {
88  r = porosity*smoothedHeaviside(epsHeaviside,phi+u) - H;
89  dr = porosity*smoothedDirac(epsDirac,phi+u);
90  }
91 
92  inline void calculateElementResidual(//element
93  double* mesh_trial_ref,
94  double* mesh_grad_trial_ref,
95  double* mesh_dof,
96  int* mesh_l2g,
97  double* dV_ref,
98  double* u_trial_ref,
99  double* u_grad_trial_ref,
100  double* u_test_ref,
101  double* u_grad_test_ref,
102  //element boundary
103  double* mesh_trial_trace_ref,
104  double* mesh_grad_trial_trace_ref,
105  double* dS_ref,
106  double* u_trial_trace_ref,
107  double* u_grad_trial_trace_ref,
108  double* u_test_trace_ref,
109  double* u_grad_test_trace_ref,
110  double* normal_ref,
111  double* boundaryJac_ref,
112  //physics
113  int nElements_global,
114  double useMetrics,
115  double epsFactHeaviside,
116  double epsFactDirac,
117  double epsFactDiffusion,
118  int* u_l2g,
119  double* elementDiameter,
120  double* nodeDiametersArray,
121  double* u_dof,
122  double* q_phi,
123  double* q_normal_phi,
124  double* ebqe_phi,
125  double* ebqe_normal_phi,
126  double* q_H,
127  double* q_u,
128  double* q_n,
129  double* ebqe_u,
130  double* ebqe_n,
131  double* q_r,
132  double* q_vos,
133  int offset_u, int stride_u,
134  double* elementResidual_u,
135  int nExteriorElementBoundaries_global,
136  int* exteriorElementBoundariesArray,
137  int* elementBoundaryElementsArray,
138  int* elementBoundaryLocalElementBoundariesArray,
139  double* element_u,
140  int eN)
141  {
142  for (int i=0;i<nDOF_test_element;i++)
143  {
144  elementResidual_u[i]=0.0;
145  }//i
146  double epsHeaviside,epsDirac,epsDiffusion,norm;
147  //loop over quadrature points and compute integrands
148  for (int k=0;k<nQuadraturePoints_element;k++)
149  {
150  //compute indeces and declare local storage
151  register int eN_k = eN*nQuadraturePoints_element+k,
152  eN_k_nSpace = eN_k*nSpace;
153  //eN_nDOF_trial_element = eN*nDOF_trial_element;
154  register double u=0.0,grad_u[nSpace],
155  r=0.0,dr=0.0,
156  jac[nSpace*nSpace],
157  jacDet,
158  jacInv[nSpace*nSpace],
159  u_grad_trial[nDOF_trial_element*nSpace],
160  u_test_dV[nDOF_trial_element],
161  u_grad_test_dV[nDOF_test_element*nSpace],
162  dV,x,y,z,
163  G[nSpace*nSpace],G_dd_G,tr_G,h_phi;
164  //
165  //compute solution and gradients at quadrature points
166  //
167  ck.calculateMapping_element(eN,
168  k,
169  mesh_dof,
170  mesh_l2g,
171  mesh_trial_ref,
172  mesh_grad_trial_ref,
173  jac,
174  jacDet,
175  jacInv,
176  x,y,z);
177  ck.calculateH_element(eN,
178  k,
179  nodeDiametersArray,
180  mesh_l2g,
181  mesh_trial_ref,
182  h_phi);
183  //get the physical integration weight
184  dV = fabs(jacDet)*dV_ref[k];
185  ck.calculateG(jacInv,G,G_dd_G,tr_G);
186 
187  /* double dir[nSpace]; */
188  /* double norm = 1.0e-8; */
189  /* for (int I=0;I<nSpace;I++) */
190  /* norm += q_normal_phi[eN_k_nSpace+I]*q_normal_phi[eN_k_nSpace+I]; */
191  /* norm = sqrt(norm); */
192  /* for (int I=0;I<nSpace;I++) */
193  /* dir[I] = q_normal_phi[eN_k_nSpace+I]/norm; */
194  /* ck.calculateGScale(G,dir,h_phi); */
195 
196  //get the trial function gradients
197  ck.gradTrialFromRef(&u_grad_trial_ref[k*nDOF_trial_element*nSpace],jacInv,u_grad_trial);
198  //get the solution
199  ck.valFromElementDOF(element_u,&u_trial_ref[k*nDOF_trial_element],u);
200  //get the solution gradients
201  ck.gradFromElementDOF(element_u,u_grad_trial,grad_u);
202  //precalculate test function products with integration weights
203  for (int j=0;j<nDOF_trial_element;j++)
204  {
205  u_test_dV[j] = u_test_ref[k*nDOF_trial_element+j]*dV;
206  for (int I=0;I<nSpace;I++)
207  {
208  u_grad_test_dV[j*nSpace+I] = u_grad_trial[j*nSpace+I]*dV;//cek warning won't work for Petrov-Galerkin
209  }
210  }
211 
212 
213 
214  //
215  //calculate pde coefficients at quadrature points
216  //
217  epsHeaviside = epsFactHeaviside*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter[eN]);
218  epsDirac = epsFactDirac* (useMetrics*h_phi+(1.0-useMetrics)*elementDiameter[eN]);
219  epsDiffusion = epsFactDiffusion*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter[eN]);
220  // *(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter[eN]);
221  evaluateCoefficients(epsHeaviside,
222  epsDirac,
223  q_phi[eN_k],
224  q_H[eN_k],
225  u,
226  1.0 - q_vos[eN_k],
227  r,
228  dr);
229  //
230  //update element residual
231  //
232  for(int i=0;i<nDOF_test_element;i++)
233  {
234  //register int eN_k_i=eN_k*nDOF_test_element+i;
235  //register int eN_k_i_nSpace = eN_k_i*nSpace;
236  register int i_nSpace=i*nSpace;
237 
238  elementResidual_u[i] += ck.Reaction_weak(r,u_test_dV[i]) +
239  ck.NumericalDiffusion(epsDiffusion,grad_u,&u_grad_test_dV[i_nSpace]);
240  }//i
241  //
242  //save momentum for time history and velocity for subgrid error
243  //save solution for other models
244  //
245 
246  q_r[eN_k] = r;
247  q_u[eN_k] = u;
248 
249 
250  norm = 1.0e-8;
251  for (int I=0;I<nSpace;I++)
252  norm += grad_u[I]*grad_u[I];
253  norm = sqrt(norm);
254  for(int I=0;I<nSpace;I++)
255  q_n[eN_k_nSpace+I] = grad_u[I]/norm;
256  }
257  }
258  void calculateResidual(arguments_dict& args)
259  {
260  xt::pyarray<double>& mesh_trial_ref = args.array<double>("mesh_trial_ref");
261  xt::pyarray<double>& mesh_grad_trial_ref = args.array<double>("mesh_grad_trial_ref");
262  xt::pyarray<double>& mesh_dof = args.array<double>("mesh_dof");
263  xt::pyarray<int>& mesh_l2g = args.array<int>("mesh_l2g");
264  xt::pyarray<double>& dV_ref = args.array<double>("dV_ref");
265  xt::pyarray<double>& u_trial_ref = args.array<double>("u_trial_ref");
266  xt::pyarray<double>& u_grad_trial_ref = args.array<double>("u_grad_trial_ref");
267  xt::pyarray<double>& u_test_ref = args.array<double>("u_test_ref");
268  xt::pyarray<double>& u_grad_test_ref = args.array<double>("u_grad_test_ref");
269  xt::pyarray<double>& mesh_trial_trace_ref = args.array<double>("mesh_trial_trace_ref");
270  xt::pyarray<double>& mesh_grad_trial_trace_ref = args.array<double>("mesh_grad_trial_trace_ref");
271  xt::pyarray<double>& dS_ref = args.array<double>("dS_ref");
272  xt::pyarray<double>& u_trial_trace_ref = args.array<double>("u_trial_trace_ref");
273  xt::pyarray<double>& u_grad_trial_trace_ref = args.array<double>("u_grad_trial_trace_ref");
274  xt::pyarray<double>& u_test_trace_ref = args.array<double>("u_test_trace_ref");
275  xt::pyarray<double>& u_grad_test_trace_ref = args.array<double>("u_grad_test_trace_ref");
276  xt::pyarray<double>& normal_ref = args.array<double>("normal_ref");
277  xt::pyarray<double>& boundaryJac_ref = args.array<double>("boundaryJac_ref");
278  int nElements_global = args.scalar<int>("nElements_global");
279  double useMetrics = args.scalar<double>("useMetrics");
280  double epsFactHeaviside = args.scalar<double>("epsFactHeaviside");
281  double epsFactDirac = args.scalar<double>("epsFactDirac");
282  double epsFactDiffusion = args.scalar<double>("epsFactDiffusion");
283  xt::pyarray<int>& u_l2g = args.array<int>("u_l2g");
284  xt::pyarray<double>& elementDiameter = args.array<double>("elementDiameter");
285  xt::pyarray<double>& nodeDiametersArray = args.array<double>("nodeDiametersArray");
286  xt::pyarray<double>& u_dof = args.array<double>("u_dof");
287  xt::pyarray<double>& q_phi = args.array<double>("q_phi");
288  xt::pyarray<double>& q_normal_phi = args.array<double>("q_normal_phi");
289  xt::pyarray<double>& ebqe_phi = args.array<double>("ebqe_phi");
290  xt::pyarray<double>& ebqe_normal_phi = args.array<double>("ebqe_normal_phi");
291  xt::pyarray<double>& q_H = args.array<double>("q_H");
292  xt::pyarray<double>& q_u = args.array<double>("q_u");
293  xt::pyarray<double>& q_n = args.array<double>("q_n");
294  xt::pyarray<double>& ebqe_u = args.array<double>("ebqe_u");
295  xt::pyarray<double>& ebqe_n = args.array<double>("ebqe_n");
296  xt::pyarray<double>& q_r = args.array<double>("q_r");
297  xt::pyarray<double>& q_vos = args.array<double>("q_vos");
298  int offset_u = args.scalar<int>("offset_u");
299  int stride_u = args.scalar<int>("stride_u");
300  xt::pyarray<double>& globalResidual = args.array<double>("globalResidual");
301  int nExteriorElementBoundaries_global = args.scalar<int>("nExteriorElementBoundaries_global");
302  xt::pyarray<int>& exteriorElementBoundariesArray = args.array<int>("exteriorElementBoundariesArray");
303  xt::pyarray<int>& elementBoundaryElementsArray = args.array<int>("elementBoundaryElementsArray");
304  xt::pyarray<int>& elementBoundaryLocalElementBoundariesArray = args.array<int>("elementBoundaryLocalElementBoundariesArray");
305  //
306  //loop over elements to compute volume integrals and load them into element and global residual
307  //
308  //eN is the element index
309  //eN_k is the quadrature point index for a scalar
310  //eN_k_nSpace is the quadrature point index for a vector
311  //eN_i is the element test function index
312  //eN_j is the element trial function index
313  //eN_k_j is the quadrature point index for a trial function
314  //eN_k_i is the quadrature point index for a trial function
315  for(int eN=0;eN<nElements_global;eN++)
316  {
317  //declare local storage for element residual and initialize
318  register double elementResidual_u[nDOF_test_element],element_u[nDOF_trial_element];
319  for (int i=0;i<nDOF_test_element;i++)
320  {
321  register int eN_i=eN*nDOF_test_element+i;
322  element_u[i] = u_dof.data()[u_l2g.data()[eN_i]];
323  }//i
324  calculateElementResidual(mesh_trial_ref.data(),
325  mesh_grad_trial_ref.data(),
326  mesh_dof.data(),
327  mesh_l2g.data(),
328  dV_ref.data(),
329  u_trial_ref.data(),
330  u_grad_trial_ref.data(),
331  u_test_ref.data(),
332  u_grad_test_ref.data(),
333  mesh_trial_trace_ref.data(),
334  mesh_grad_trial_trace_ref.data(),
335  dS_ref.data(),
336  u_trial_trace_ref.data(),
337  u_grad_trial_trace_ref.data(),
338  u_test_trace_ref.data(),
339  u_grad_test_trace_ref.data(),
340  normal_ref.data(),
341  boundaryJac_ref.data(),
342  nElements_global,
343  useMetrics,
344  epsFactHeaviside,
345  epsFactDirac,
346  epsFactDiffusion,
347  u_l2g.data(),
348  elementDiameter.data(),
349  nodeDiametersArray.data(),
350  u_dof.data(),
351  q_phi.data(),
352  q_normal_phi.data(),
353  ebqe_phi.data(),
354  ebqe_normal_phi.data(),
355  q_H.data(),
356  q_u.data(),
357  q_n.data(),
358  ebqe_u.data(),
359  ebqe_n.data(),
360  q_r.data(),
361  q_vos.data(),
362  offset_u,stride_u,
363  elementResidual_u,
364  nExteriorElementBoundaries_global,
365  exteriorElementBoundariesArray.data(),
366  elementBoundaryElementsArray.data(),
367  elementBoundaryLocalElementBoundariesArray.data(),
368  element_u,
369  eN);
370  //
371  //load element into global residual and save element residual
372  //
373  for(int i=0;i<nDOF_test_element;i++)
374  {
375  register int eN_i=eN*nDOF_test_element+i;
376 
377  globalResidual.data()[offset_u+stride_u*u_l2g.data()[eN_i]]+=elementResidual_u[i];
378  }//i
379  }//elements
380  //
381  //loop over exterior element boundaries to calculate levelset gradient
382  //
383  //ebNE is the Exterior element boundary INdex
384  //ebN is the element boundary INdex
385  //eN is the element index
386  for (int ebNE = 0; ebNE < nExteriorElementBoundaries_global; ebNE++)
387  {
388  register int ebN = exteriorElementBoundariesArray.data()[ebNE],
389  eN = elementBoundaryElementsArray.data()[ebN*2+0],
390  ebN_local = elementBoundaryLocalElementBoundariesArray.data()[ebN*2+0];
391  //eN_nDOF_trial_element = eN*nDOF_trial_element;
392  //register double elementResidual_u[nDOF_test_element];
393  double element_u[nDOF_trial_element];
394  for (int i=0;i<nDOF_test_element;i++)
395  {
396  register int eN_i=eN*nDOF_test_element+i;
397  element_u[i] = u_dof.data()[u_l2g.data()[eN_i]];
398  }//i
399  for (int kb=0;kb<nQuadraturePoints_elementBoundary;kb++)
400  {
401  register int ebNE_kb = ebNE*nQuadraturePoints_elementBoundary+kb,
402  ebNE_kb_nSpace = ebNE_kb*nSpace,
403  ebN_local_kb = ebN_local*nQuadraturePoints_elementBoundary+kb,
404  ebN_local_kb_nSpace = ebN_local_kb*nSpace;
405  register double u_ext=0.0,
406  grad_u_ext[nSpace],
407  jac_ext[nSpace*nSpace],
408  jacDet_ext,
409  jacInv_ext[nSpace*nSpace],
410  boundaryJac[nSpace*(nSpace-1)],
411  metricTensor[(nSpace-1)*(nSpace-1)],
412  metricTensorDetSqrt,
413  u_grad_trial_trace[nDOF_trial_element*nSpace],
414  normal[nSpace],x_ext,y_ext,z_ext,
415  G[nSpace*nSpace],G_dd_G,tr_G,norm;
416  //
417  //calculate the solution and gradients at quadrature points
418  //
419  ck.calculateMapping_elementBoundary(eN,
420  ebN_local,
421  kb,
422  ebN_local_kb,
423  mesh_dof.data(),
424  mesh_l2g.data(),
425  mesh_trial_trace_ref.data(),
426  mesh_grad_trial_trace_ref.data(),
427  boundaryJac_ref.data(),
428  jac_ext,
429  jacDet_ext,
430  jacInv_ext,
431  boundaryJac,
432  metricTensor,
433  metricTensorDetSqrt,
434  normal_ref.data(),
435  normal,
436  x_ext,y_ext,z_ext);
437  //get the metric tensor
438  //cek todo use symmetry
439  ck.calculateG(jacInv_ext,G,G_dd_G,tr_G);
440  //compute shape and solution information
441  //shape
442  ck.gradTrialFromRef(&u_grad_trial_trace_ref.data()[ebN_local_kb_nSpace*nDOF_trial_element],jacInv_ext,u_grad_trial_trace);
443  //solution and gradients
444  ck.valFromElementDOF(element_u,&u_trial_trace_ref.data()[ebN_local_kb*nDOF_test_element],u_ext);
445  ck.gradFromElementDOF(element_u,u_grad_trial_trace,grad_u_ext);
446 
447  ebqe_u.data()[ebNE_kb] = u_ext;
448  norm = 1.0e-8;
449  for (int I=0;I<nSpace;I++)
450  norm += grad_u_ext[I]*grad_u_ext[I];
451  norm = sqrt(norm);
452  for (int I=0;I<nSpace;I++)
453  ebqe_n.data()[ebNE_kb_nSpace+I] = grad_u_ext[I]/norm;
454  }//kb
455  }//ebNE
456 
457 
458  }
459 
460  inline void calculateElementJacobian(//element
461  double* mesh_trial_ref,
462  double* mesh_grad_trial_ref,
463  double* mesh_dof,
464  int* mesh_l2g,
465  double* dV_ref,
466  double* u_trial_ref,
467  double* u_grad_trial_ref,
468  double* u_test_ref,
469  double* u_grad_test_ref,
470  //element boundary
471  double* mesh_trial_trace_ref,
472  double* mesh_grad_trial_trace_ref,
473  double* dS_ref,
474  double* u_trial_trace_ref,
475  double* u_grad_trial_trace_ref,
476  double* u_test_trace_ref,
477  double* u_grad_test_trace_ref,
478  double* normal_ref,
479  double* boundaryJac_ref,
480  //physics
481  int nElements_global,
482  double useMetrics,
483  double epsFactHeaviside,
484  double epsFactDirac,
485  double epsFactDiffusion,
486  int* u_l2g,
487  double* elementDiameter,
488  double* nodeDiametersArray,
489  double* u_dof,
490  // double* u_trial,
491  // double* u_grad_trial,
492  // double* u_test_dV,
493  // double* u_grad_test_dV,
494  double* q_phi,
495  double* q_normal_phi,
496  double* q_H,
497  double* q_vos,
498  double* elementJacobian_u_u,
499  double* element_u,
500  int eN)
501  {
502  for (int i=0;i<nDOF_test_element;i++)
503  for (int j=0;j<nDOF_trial_element;j++)
504  {
505  elementJacobian_u_u[i*nDOF_trial_element+j]=0.0;
506  }
507  double epsHeaviside,epsDirac,epsDiffusion;
508  for (int k=0;k<nQuadraturePoints_element;k++)
509  {
510  int eN_k = eN*nQuadraturePoints_element+k; //index to a scalar at a quadrature point
511 
512  //declare local storage
513  register double u=0.0,
514  grad_u[nSpace],
515  r=0.0,dr=0.0,
516  jac[nSpace*nSpace],
517  jacDet,
518  jacInv[nSpace*nSpace],
519  u_grad_trial[nDOF_trial_element*nSpace],
520  dV,
521  u_test_dV[nDOF_test_element],
522  u_grad_test_dV[nDOF_test_element*nSpace],
523  x,y,z,
524  G[nSpace*nSpace],G_dd_G,tr_G,h_phi;
525  //
526  //calculate solution and gradients at quadrature points
527  //
528  ck.calculateMapping_element(eN,
529  k,
530  mesh_dof,
531  mesh_l2g,
532  mesh_trial_ref,
533  mesh_grad_trial_ref,
534  jac,
535  jacDet,
536  jacInv,
537  x,y,z);
538  ck.calculateH_element(eN,
539  k,
540  nodeDiametersArray,
541  mesh_l2g,
542  mesh_trial_ref,
543  h_phi);
544  //get the physical integration weight
545  dV = fabs(jacDet)*dV_ref[k];
546  ck.calculateG(jacInv,G,G_dd_G,tr_G);
547 
548  /* double dir[nSpace]; */
549  /* double norm = 1.0e-8; */
550  /* for (int I=0;I<nSpace;I++) */
551  /* norm += q_normal_phi[eN_k_nSpace+I]*q_normal_phi[eN_k_nSpace+I]; */
552  /* norm = sqrt(norm); */
553  /* for (int I=0;I<nSpace;I++) */
554  /* dir[I] = q_normal_phi[eN_k_nSpace+I]/norm; */
555  /* ck.calculateGScale(G,dir,h_phi); */
556 
557 
558  //get the trial function gradients
559  ck.gradTrialFromRef(&u_grad_trial_ref[k*nDOF_trial_element*nSpace],jacInv,u_grad_trial);
560  //get the solution
561  ck.valFromElementDOF(element_u,&u_trial_ref[k*nDOF_trial_element],u);
562  //get the solution gradients
563  ck.gradFromElementDOF(element_u,u_grad_trial,grad_u);
564  //precalculate test function products with integration weights
565  for (int j=0;j<nDOF_trial_element;j++)
566  {
567  u_test_dV[j] = u_test_ref[k*nDOF_trial_element+j]*dV;
568  for (int I=0;I<nSpace;I++)
569  {
570  u_grad_test_dV[j*nSpace+I] = u_grad_trial[j*nSpace+I]*dV;//cek warning won't work for Petrov-Galerkin
571  }
572  }
573  //
574  //calculate pde coefficients and derivatives at quadrature points
575  //
576  epsHeaviside=epsFactHeaviside*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter[eN]);
577  epsDirac =epsFactDirac* (useMetrics*h_phi+(1.0-useMetrics)*elementDiameter[eN]);
578  epsDiffusion=epsFactDiffusion*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter[eN]);
579  // *(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter[eN]);
580  evaluateCoefficients(epsHeaviside,
581  epsDirac,
582  q_phi[eN_k],
583  q_H[eN_k],
584  u,
585  1.0 - q_vos[eN_k],
586  r,
587  dr);
588  for(int i=0;i<nDOF_test_element;i++)
589  {
590  //int eN_k_i=eN_k*nDOF_test_element+i;
591  //int eN_k_i_nSpace=eN_k_i*nSpace;
592  int i_nSpace=i*nSpace;
593  for(int j=0;j<nDOF_trial_element;j++)
594  {
595  //int eN_k_j=eN_k*nDOF_trial_element+j;
596  //int eN_k_j_nSpace = eN_k_j*nSpace;
597  int j_nSpace = j*nSpace;
598 
599  elementJacobian_u_u[i*nDOF_trial_element+j] += ck.ReactionJacobian_weak(dr,u_trial_ref[k*nDOF_trial_element+j],u_test_dV[i]) +
600  ck.NumericalDiffusionJacobian(epsDiffusion,&u_grad_trial[j_nSpace],&u_grad_test_dV[i_nSpace]);
601  }//j
602  }//i
603  }//k
604  }
605  void calculateJacobian(arguments_dict& args)
606  {
607  xt::pyarray<double>& mesh_trial_ref = args.array<double>("mesh_trial_ref");
608  xt::pyarray<double>& mesh_grad_trial_ref = args.array<double>("mesh_grad_trial_ref");
609  xt::pyarray<double>& mesh_dof = args.array<double>("mesh_dof");
610  xt::pyarray<int>& mesh_l2g = args.array<int>("mesh_l2g");
611  xt::pyarray<double>& dV_ref = args.array<double>("dV_ref");
612  xt::pyarray<double>& u_trial_ref = args.array<double>("u_trial_ref");
613  xt::pyarray<double>& u_grad_trial_ref = args.array<double>("u_grad_trial_ref");
614  xt::pyarray<double>& u_test_ref = args.array<double>("u_test_ref");
615  xt::pyarray<double>& u_grad_test_ref = args.array<double>("u_grad_test_ref");
616  xt::pyarray<double>& mesh_trial_trace_ref = args.array<double>("mesh_trial_trace_ref");
617  xt::pyarray<double>& mesh_grad_trial_trace_ref = args.array<double>("mesh_grad_trial_trace_ref");
618  xt::pyarray<double>& dS_ref = args.array<double>("dS_ref");
619  xt::pyarray<double>& u_trial_trace_ref = args.array<double>("u_trial_trace_ref");
620  xt::pyarray<double>& u_grad_trial_trace_ref = args.array<double>("u_grad_trial_trace_ref");
621  xt::pyarray<double>& u_test_trace_ref = args.array<double>("u_test_trace_ref");
622  xt::pyarray<double>& u_grad_test_trace_ref = args.array<double>("u_grad_test_trace_ref");
623  xt::pyarray<double>& normal_ref = args.array<double>("normal_ref");
624  xt::pyarray<double>& boundaryJac_ref = args.array<double>("boundaryJac_ref");
625  int nElements_global = args.scalar<int>("nElements_global");
626  double useMetrics = args.scalar<double>("useMetrics");
627  double epsFactHeaviside = args.scalar<double>("epsFactHeaviside");
628  double epsFactDirac = args.scalar<double>("epsFactDirac");
629  double epsFactDiffusion = args.scalar<double>("epsFactDiffusion");
630  xt::pyarray<int>& u_l2g = args.array<int>("u_l2g");
631  xt::pyarray<double>& elementDiameter = args.array<double>("elementDiameter");
632  xt::pyarray<double>& nodeDiametersArray = args.array<double>("nodeDiametersArray");
633  xt::pyarray<double>& u_dof = args.array<double>("u_dof");
634  xt::pyarray<double>& q_phi = args.array<double>("q_phi");
635  xt::pyarray<double>& q_normal_phi = args.array<double>("q_normal_phi");
636  xt::pyarray<double>& q_H = args.array<double>("q_H");
637  xt::pyarray<double>& q_vos = args.array<double>("q_vos");
638  xt::pyarray<int>& csrRowIndeces_u_u = args.array<int>("csrRowIndeces_u_u");
639  xt::pyarray<int>& csrColumnOffsets_u_u = args.array<int>("csrColumnOffsets_u_u");
640  xt::pyarray<double>& globalJacobian = args.array<double>("globalJacobian");
641  //
642  //loop over elements to compute volume integrals and load them into the element Jacobians and global Jacobian
643  //
644  for(int eN=0;eN<nElements_global;eN++)
645  {
646  register double elementJacobian_u_u[nDOF_test_element*nDOF_trial_element],element_u[nDOF_trial_element];
647  for (int j=0;j<nDOF_trial_element;j++)
648  {
649  register int eN_j = eN*nDOF_trial_element+j;
650  element_u[j] = u_dof.data()[u_l2g.data()[eN_j]];
651  }
652  calculateElementJacobian(mesh_trial_ref.data(),
653  mesh_grad_trial_ref.data(),
654  mesh_dof.data(),
655  mesh_l2g.data(),
656  dV_ref.data(),
657  u_trial_ref.data(),
658  u_grad_trial_ref.data(),
659  u_test_ref.data(),
660  u_grad_test_ref.data(),
661  mesh_trial_trace_ref.data(),
662  mesh_grad_trial_trace_ref.data(),
663  dS_ref.data(),
664  u_trial_trace_ref.data(),
665  u_grad_trial_trace_ref.data(),
666  u_test_trace_ref.data(),
667  u_grad_test_trace_ref.data(),
668  normal_ref.data(),
669  boundaryJac_ref.data(),
670  nElements_global,
671  useMetrics,
672  epsFactHeaviside,
673  epsFactDirac,
674  epsFactDiffusion,
675  u_l2g.data(),
676  elementDiameter.data(),
677  nodeDiametersArray.data(),
678  u_dof.data(),
679  q_phi.data(),
680  q_normal_phi.data(),
681  q_H.data(),
682  q_vos.data(),
683  elementJacobian_u_u,
684  element_u,
685  eN);
686  //
687  //load into element Jacobian into global Jacobian
688  //
689  for (int i=0;i<nDOF_test_element;i++)
690  {
691  int eN_i = eN*nDOF_test_element+i;
692  for (int j=0;j<nDOF_trial_element;j++)
693  {
694  int eN_i_j = eN_i*nDOF_trial_element+j;
695 
696  globalJacobian.data()[csrRowIndeces_u_u[eN_i] + csrColumnOffsets_u_u[eN_i_j]] += elementJacobian_u_u[i*nDOF_trial_element+j];
697  }//j
698  }//i
699  }//elements
700  }//computeJacobian
701  };//cppPresInit
702 
703  inline cppPresInit_base* newPresInit(int nSpaceIn,
704  int nQuadraturePoints_elementIn,
705  int nDOF_mesh_trial_elementIn,
706  int nDOF_trial_elementIn,
707  int nDOF_test_elementIn,
708  int nQuadraturePoints_elementBoundaryIn,
709  int CompKernelFlag)
710  {
711  if (nSpaceIn == 2)
712  return proteus::chooseAndAllocateDiscretization2D<cppPresInit_base,cppPresInit,CompKernel>(nSpaceIn,
713  nQuadraturePoints_elementIn,
714  nDOF_mesh_trial_elementIn,
715  nDOF_trial_elementIn,
716  nDOF_test_elementIn,
717  nQuadraturePoints_elementBoundaryIn,
718  CompKernelFlag);
719  else
720  return proteus::chooseAndAllocateDiscretization<cppPresInit_base,cppPresInit,CompKernel>(nSpaceIn,
721  nQuadraturePoints_elementIn,
722  nDOF_mesh_trial_elementIn,
723  nDOF_trial_elementIn,
724  nDOF_test_elementIn,
725  nQuadraturePoints_elementBoundaryIn,
726  CompKernelFlag);
727  }
728 }//proteus
729 #endif
proteus::cppPresInit::calculateResidual
void calculateResidual(arguments_dict &args)
Definition: PresInit.h:258
proteus::cppPresInit_base::calculateResidual
virtual void calculateResidual(arguments_dict &args)=0
HI
#define HI
Definition: Headers.h:4
proteus::cppPresInit::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, double useMetrics, double epsFactHeaviside, double epsFactDirac, double epsFactDiffusion, int *u_l2g, double *elementDiameter, double *nodeDiametersArray, double *u_dof, double *q_phi, double *q_normal_phi, double *ebqe_phi, double *ebqe_normal_phi, double *q_H, double *q_u, double *q_n, double *ebqe_u, double *ebqe_n, double *q_r, double *q_vos, int offset_u, int stride_u, double *elementResidual_u, int nExteriorElementBoundaries_global, int *exteriorElementBoundariesArray, int *elementBoundaryElementsArray, int *elementBoundaryLocalElementBoundariesArray, double *element_u, int eN)
Definition: PresInit.h:92
proteus::cppPresInit_base
Definition: PresInit.h:15
proteus::cppPresInit::smoothedDirac
double smoothedDirac(double eps, double phi)
Definition: PresInit.h:67
proteus::cppPresInit::smoothedHeaviside
double smoothedHeaviside(double eps, double phi)
Definition: PresInit.h:35
ModelFactory.h
CompKernel.h
proteus::arguments_dict::scalar
T & scalar(const std::string &key)
proteus::cppPresInit::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, double useMetrics, double epsFactHeaviside, double epsFactDirac, double epsFactDiffusion, int *u_l2g, double *elementDiameter, double *nodeDiametersArray, double *u_dof, double *q_phi, double *q_normal_phi, double *q_H, double *q_vos, double *elementJacobian_u_u, double *element_u, int eN)
Definition: PresInit.h:460
proteus::arguments_dict::array
xt::pyarray< T > & array(const std::string &key)
H
Double H
Definition: Headers.h:65
proteus::cppPresInit::calculateJacobian
void calculateJacobian(arguments_dict &args)
Definition: PresInit.h:605
proteus::cppPresInit::evaluateCoefficients
void evaluateCoefficients(const double &epsHeaviside, const double &epsDirac, const double &phi, const double &H, const double &u, const double &porosity, double &r, double &dr)
Definition: PresInit.h:79
proteus::cppPresInit_base::calculateJacobian
virtual void calculateJacobian(arguments_dict &args)=0
z
Double * z
Definition: Headers.h:49
u
Double u
Definition: Headers.h:89
proteus::phi
double phi(const double &g, const double &h, const double &hL, const double &hR, const double &uL, const double &uR)
Definition: SW2DCV.h:62
proteus::cppPresInit_base::~cppPresInit_base
virtual ~cppPresInit_base()
Definition: PresInit.h:17
proteus
Definition: ADR.h:17
proteus::cppPresInit
Definition: PresInit.h:30
proteus::cppPresInit::ck
CompKernelType ck
Definition: PresInit.h:32
proteus::newPresInit
cppPresInit_base * newPresInit(int nSpaceIn, int nQuadraturePoints_elementIn, int nDOF_mesh_trial_elementIn, int nDOF_trial_elementIn, int nDOF_test_elementIn, int nQuadraturePoints_elementBoundaryIn, int CompKernelFlag)
Definition: PresInit.h:703
r
Double r
Definition: Headers.h:83
proteus::arguments_dict
Definition: ArgumentsDict.h:70
proteus::cppPresInit::smoothedHeaviside_integral
double smoothedHeaviside_integral(double eps, double phi)
Definition: PresInit.h:49
proteus::cppPresInit::cppPresInit
cppPresInit()
Definition: PresInit.h:33
ArgumentsDict.h
Generated on Fri Jul 1 2022 11:17:39 for proteus by doxygen 1.8.20