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