proteus  1.8.1
C/C++/Fortran libraries
MCorr.h
Go to the documentation of this file.
1 #ifndef MCorr_H
2 #define MCorr_H
3 #include <cmath>
4 #include <iostream>
5 #include <valarray>
6 #include "CompKernel.h"
7 #include "ModelFactory.h"
8 #include "equivalent_polynomials.h"
9 #include PROTEUS_LAPACK_H
10 #include "ArgumentsDict.h"
11 #include "xtensor-python/pyarray.hpp"
12 
13 namespace py = pybind11;
14 
15 namespace proteus
16 {
17 
18  template<int nSpace, int nP, int nQ, int nEBQ>
19  using GeneralizedFunctions = equivalent_polynomials::GeneralizedFunctions_mix<nSpace, nP, nQ, nEBQ>;
20  //using GeneralizedFunctions = equivalent_polynomials::Regularized<nSpace, nP, nQ>;
21  //using GeneralizedFunctions = equivalent_polynomials::EquivalentPolynomials<nSpace, nP, nQ>;
22 
23  class MCorr_base
24  {
25  public:
26  std::valarray<double> Rpos, Rneg;
27  std::valarray<double> FluxCorrectionMatrix;
28  virtual ~MCorr_base(){}
29  virtual void calculateResidual(arguments_dict& args, bool useExact)=0;
30  virtual void calculateJacobian(arguments_dict& args, bool useExact)=0;
31  virtual void elementSolve(arguments_dict& args)=0;
32  virtual void elementConstantSolve(arguments_dict& args)=0;
33  virtual std::tuple<double, double> globalConstantRJ(arguments_dict& args)=0;
34  virtual double calculateMass(arguments_dict& args, bool useExact)=0;
35  virtual void setMassQuadrature(arguments_dict& args, bool useExact)=0;
36  virtual void FCTStep(arguments_dict& args)=0;
37  virtual void calculateMassMatrix(arguments_dict& args)=0;
38  virtual void setMassQuadratureEdgeBasedStabilizationMethods(arguments_dict& args,
39  bool useExact)=0;
40  };
41 
42  template<class CompKernelType,
43  int nSpace,
44  int nQuadraturePoints_element,
45  int nDOF_mesh_trial_element,
46  int nDOF_trial_element,
47  int nDOF_test_element,
48  int nQuadraturePoints_elementBoundary>
49  class MCorr : public MCorr_base
50  {
51  public:
52  CompKernelType ck;
53  GeneralizedFunctions<nSpace,2,nQuadraturePoints_element,nQuadraturePoints_elementBoundary> gf;
54  GeneralizedFunctions<nSpace,2,nDOF_trial_element,nQuadraturePoints_elementBoundary> gf_nodes;
55  MCorr():ck()
56  {}
57 
58  inline
59  void evaluateCoefficients(const double& epsHeaviside,
60  const double& epsDirac,
61  const double& phi,
62  const double& H,
63  const double& u,
64  const double& porosity,
65  double& r,
66  double& dr)
67  {
68  r = porosity*(gf.H(epsHeaviside,phi+u) - H);
69  dr = porosity*gf.D(epsDirac,phi+u);
70  }
71 
72  inline void calculateElementResidual(//element
73  double* mesh_trial_ref,
74  double* mesh_grad_trial_ref,
75  double* mesh_dof,
76  int* mesh_l2g,
77  double* dV_ref,
78  double* u_trial_ref,
79  double* u_grad_trial_ref,
80  double* u_test_ref,
81  double* u_grad_test_ref,
82  //element boundary
83  double* mesh_trial_trace_ref,
84  double* mesh_grad_trial_trace_ref,
85  double* dS_ref,
86  double* u_trial_trace_ref,
87  double* u_grad_trial_trace_ref,
88  double* u_test_trace_ref,
89  double* u_grad_test_trace_ref,
90  double* normal_ref,
91  double* boundaryJac_ref,
92  //physics
93  int nElements_global,
94  double useMetrics,
95  double epsFactHeaviside,
96  double epsFactDirac,
97  double epsFactDiffusion,
98  int* u_l2g,
99  double* elementDiameter,
100  double* nodeDiametersArray,
101  double* u_dof,
102  double* q_phi,
103  double* q_normal_phi,
104  double* ebqe_phi,
105  double* ebqe_normal_phi,
106  double* q_H,
107  double* q_u,
108  double* q_n,
109  double* ebqe_u,
110  double* ebqe_n,
111  double* q_r,
112  double* q_porosity,
113  int offset_u, int stride_u,
114  double* elementResidual_u,
115  int nExteriorElementBoundaries_global,
116  int* exteriorElementBoundariesArray,
117  int* elementBoundaryElementsArray,
118  int* elementBoundaryLocalElementBoundariesArray,
119  double* element_u,
120  int eN)
121  {
122  for (int i=0;i<nDOF_test_element;i++)
123  {
124  elementResidual_u[i]=0.0;
125  }//i
126  double epsHeaviside,epsDirac,epsDiffusion,norm;
127  //loop over quadrature points and compute integrands
128  for (int k=0;k<nQuadraturePoints_element;k++)
129  {
130  //compute indeces and declare local storage
131  register int eN_k = eN*nQuadraturePoints_element+k,
132  eN_k_nSpace = eN_k*nSpace;
133  //eN_nDOF_trial_element = eN*nDOF_trial_element;
134  register double u=0.0,grad_u[nSpace],
135  r=0.0,dr=0.0,
136  jac[nSpace*nSpace],
137  jacDet,
138  jacInv[nSpace*nSpace],
139  u_grad_trial[nDOF_trial_element*nSpace],
140  u_test_dV[nDOF_trial_element],
141  u_grad_test_dV[nDOF_test_element*nSpace],
142  dV,x,y,z,
143  G[nSpace*nSpace],G_dd_G,tr_G,h_phi;
144  gf.set_quad(k);
145  //
146  //compute solution and gradients at quadrature points
147  //
148  ck.calculateMapping_element(eN,
149  k,
150  mesh_dof,
151  mesh_l2g,
152  mesh_trial_ref,
153  mesh_grad_trial_ref,
154  jac,
155  jacDet,
156  jacInv,
157  x,y,z);
158  ck.calculateH_element(eN,
159  k,
160  nodeDiametersArray,
161  mesh_l2g,
162  mesh_trial_ref,
163  h_phi);
164  //get the physical integration weight
165  dV = fabs(jacDet)*dV_ref[k];
166  ck.calculateG(jacInv,G,G_dd_G,tr_G);
167 
168  /* double dir[nSpace]; */
169  /* double norm = 1.0e-8; */
170  /* for (int I=0;I<nSpace;I++) */
171  /* norm += q_normal_phi[eN_k_nSpace+I]*q_normal_phi[eN_k_nSpace+I]; */
172  /* norm = sqrt(norm); */
173  /* for (int I=0;I<nSpace;I++) */
174  /* dir[I] = q_normal_phi[eN_k_nSpace+I]/norm; */
175  /* ck.calculateGScale(G,dir,h_phi); */
176 
177  //get the trial function gradients
178  ck.gradTrialFromRef(&u_grad_trial_ref[k*nDOF_trial_element*nSpace],jacInv,u_grad_trial);
179  //get the solution
180  ck.valFromElementDOF(element_u,&u_trial_ref[k*nDOF_trial_element],u);
181  //get the solution gradients
182  ck.gradFromElementDOF(element_u,u_grad_trial,grad_u);
183  //precalculate test function products with integration weights
184  for (int j=0;j<nDOF_trial_element;j++)
185  {
186  u_test_dV[j] = u_test_ref[k*nDOF_trial_element+j]*dV;
187  for (int I=0;I<nSpace;I++)
188  {
189  u_grad_test_dV[j*nSpace+I] = u_grad_trial[j*nSpace+I]*dV;//cek warning won't work for Petrov-Galerkin
190  }
191  }
192 
193 
194 
195  //
196  //calculate pde coefficients at quadrature points
197  //
198  epsHeaviside = epsFactHeaviside*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter[eN]);
199  epsDirac = epsFactDirac* (useMetrics*h_phi+(1.0-useMetrics)*elementDiameter[eN]);
200  epsDiffusion = epsFactDiffusion*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter[eN]);
201  // *(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter[eN]);
202  evaluateCoefficients(epsHeaviside,
203  epsDirac,
204  q_phi[eN_k],
205  q_H[eN_k],
206  u,
207  q_porosity[eN_k],
208  r,
209  dr);
210  //
211  //update element residual
212  //
213  for(int i=0;i<nDOF_test_element;i++)
214  {
215  //register int eN_k_i=eN_k*nDOF_test_element+i;
216  //register int eN_k_i_nSpace = eN_k_i*nSpace;
217  register int i_nSpace=i*nSpace;
218 
219  elementResidual_u[i] += ck.Reaction_weak(r,u_test_dV[i]) +
220  ck.NumericalDiffusion(epsDiffusion,grad_u,&u_grad_test_dV[i_nSpace]);
221  }//i
222  //
223  //save momentum for time history and velocity for subgrid error
224  //save solution for other models
225  //
226 
227  q_r[eN_k] = r;
228  q_u[eN_k] = u;
229 
230 
231  norm = 1.0e-8;
232  for (int I=0;I<nSpace;I++)
233  norm += grad_u[I]*grad_u[I];
234  norm = sqrt(norm);
235  for(int I=0;I<nSpace;I++)
236  q_n[eN_k_nSpace+I] = grad_u[I]/norm;
237  }
238  }
239  void calculateResidual(arguments_dict& args,
240  bool useExact)
241  {
242  xt::pyarray<double>& mesh_trial_ref = args.array<double>("mesh_trial_ref");
243  xt::pyarray<double>& mesh_grad_trial_ref = args.array<double>("mesh_grad_trial_ref");
244  xt::pyarray<double>& mesh_dof = args.array<double>("mesh_dof");
245  xt::pyarray<int>& mesh_l2g = args.array<int>("mesh_l2g");
246  xt::pyarray<double>& x_ref = args.array<double>("x_ref");
247  xt::pyarray<double>& dV_ref = args.array<double>("dV_ref");
248  xt::pyarray<double>& u_trial_ref = args.array<double>("u_trial_ref");
249  xt::pyarray<double>& u_grad_trial_ref = args.array<double>("u_grad_trial_ref");
250  xt::pyarray<double>& u_test_ref = args.array<double>("u_test_ref");
251  xt::pyarray<double>& u_grad_test_ref = args.array<double>("u_grad_test_ref");
252  xt::pyarray<double>& mesh_trial_trace_ref = args.array<double>("mesh_trial_trace_ref");
253  xt::pyarray<double>& mesh_grad_trial_trace_ref = args.array<double>("mesh_grad_trial_trace_ref");
254  xt::pyarray<double>& dS_ref = args.array<double>("dS_ref");
255  xt::pyarray<double>& u_trial_trace_ref = args.array<double>("u_trial_trace_ref");
256  xt::pyarray<double>& u_grad_trial_trace_ref = args.array<double>("u_grad_trial_trace_ref");
257  xt::pyarray<double>& u_test_trace_ref = args.array<double>("u_test_trace_ref");
258  xt::pyarray<double>& u_grad_test_trace_ref = args.array<double>("u_grad_test_trace_ref");
259  xt::pyarray<double>& normal_ref = args.array<double>("normal_ref");
260  xt::pyarray<double>& boundaryJac_ref = args.array<double>("boundaryJac_ref");
261  int nElements_global = args.scalar<int>("nElements_global");
262  double useMetrics = args.scalar<double>("useMetrics");
263  double epsFactHeaviside = args.scalar<double>("epsFactHeaviside");
264  double epsFactDirac = args.scalar<double>("epsFactDirac");
265  double epsFactDiffusion = args.scalar<double>("epsFactDiffusion");
266  xt::pyarray<int>& u_l2g = args.array<int>("u_l2g");
267  xt::pyarray<int>& r_l2g = args.array<int>("r_l2g");
268  xt::pyarray<double>& elementDiameter = args.array<double>("elementDiameter");
269  xt::pyarray<double>& nodeDiametersArray = args.array<double>("nodeDiametersArray");
270  xt::pyarray<double>& u_dof = args.array<double>("u_dof");
271  xt::pyarray<double>& phi_dof = args.array<double>("phi_dof");
272  xt::pyarray<double>& q_phi = args.array<double>("q_phi");
273  xt::pyarray<double>& q_normal_phi = args.array<double>("q_normal_phi");
274  xt::pyarray<double>& ebqe_phi = args.array<double>("ebqe_phi");
275  xt::pyarray<double>& ebqe_normal_phi = args.array<double>("ebqe_normal_phi");
276  xt::pyarray<double>& q_H = args.array<double>("q_H");
277  xt::pyarray<double>& q_u = args.array<double>("q_u");
278  xt::pyarray<double>& q_n = args.array<double>("q_n");
279  xt::pyarray<double>& ebqe_u = args.array<double>("ebqe_u");
280  xt::pyarray<double>& ebqe_n = args.array<double>("ebqe_n");
281  xt::pyarray<double>& q_r = args.array<double>("q_r");
282  xt::pyarray<double>& q_porosity = args.array<double>("q_porosity");
283  int offset_u = args.scalar<int>("offset_u");
284  int stride_u = args.scalar<int>("stride_u");
285  xt::pyarray<double>& globalResidual = args.array<double>("globalResidual");
286  int nExteriorElementBoundaries_global = args.scalar<int>("nExteriorElementBoundaries_global");
287  xt::pyarray<int>& exteriorElementBoundariesArray = args.array<int>("exteriorElementBoundariesArray");
288  xt::pyarray<int>& elementBoundaryElementsArray = args.array<int>("elementBoundaryElementsArray");
289  xt::pyarray<int>& elementBoundaryLocalElementBoundariesArray = args.array<int>("elementBoundaryLocalElementBoundariesArray");
290  //
291  //loop over elements to compute volume integrals and load them into element and global residual
292  //
293  //eN is the element index
294  //eN_k is the quadrature point index for a scalar
295  //eN_k_nSpace is the quadrature point index for a vector
296  //eN_i is the element test function index
297  //eN_j is the element trial function index
298  //eN_k_j is the quadrature point index for a trial function
299  //eN_k_i is the quadrature point index for a trial function
300  gf.useExact = useExact;
301  for(int eN=0;eN<nElements_global;eN++)
302  {
303  //declare local storage for element residual and initialize
304  register double elementResidual_u[nDOF_test_element],element_u[nDOF_trial_element],element_phi[nDOF_trial_element];
305  for (int i=0;i<nDOF_test_element;i++)
306  {
307  register int eN_i=eN*nDOF_test_element+i;
308  element_u[i] = u_dof.data()[u_l2g.data()[eN_i]];
309  element_phi[i] = phi_dof.data()[u_l2g.data()[eN_i]] + element_u[i];
310  }//i
311  double element_nodes[nDOF_mesh_trial_element*3];
312  for (int i=0;i<nDOF_mesh_trial_element;i++)
313  {
314  register int eN_i=eN*nDOF_mesh_trial_element+i;
315  for(int I=0;I<3;I++)
316  element_nodes[i*3 + I] = mesh_dof.data()[mesh_l2g.data()[eN_i]*3 + I];
317  }//i
318  gf.calculate(element_phi, element_nodes, x_ref.data(),false);
319  calculateElementResidual(mesh_trial_ref.data(),
320  mesh_grad_trial_ref.data(),
321  mesh_dof.data(),
322  mesh_l2g.data(),
323  dV_ref.data(),
324  u_trial_ref.data(),
325  u_grad_trial_ref.data(),
326  u_test_ref.data(),
327  u_grad_test_ref.data(),
328  mesh_trial_trace_ref.data(),
329  mesh_grad_trial_trace_ref.data(),
330  dS_ref.data(),
331  u_trial_trace_ref.data(),
332  u_grad_trial_trace_ref.data(),
333  u_test_trace_ref.data(),
334  u_grad_test_trace_ref.data(),
335  normal_ref.data(),
336  boundaryJac_ref.data(),
337  nElements_global,
338  useMetrics,
339  epsFactHeaviside,
340  epsFactDirac,
341  epsFactDiffusion,
342  u_l2g.data(),
343  elementDiameter.data(),
344  nodeDiametersArray.data(),
345  u_dof.data(),
346  q_phi.data(),
347  q_normal_phi.data(),
348  ebqe_phi.data(),
349  ebqe_normal_phi.data(),
350  q_H.data(),
351  q_u.data(),
352  q_n.data(),
353  ebqe_u.data(),
354  ebqe_n.data(),
355  q_r.data(),
356  q_porosity.data(),
357  offset_u,stride_u,
358  elementResidual_u,
359  nExteriorElementBoundaries_global,
360  exteriorElementBoundariesArray.data(),
361  elementBoundaryElementsArray.data(),
362  elementBoundaryLocalElementBoundariesArray.data(),
363  element_u,
364  eN);
365  //
366  //load element into global residual and save element residual
367  //
368  for(int i=0;i<nDOF_test_element;i++)
369  {
370  register int eN_i=eN*nDOF_test_element+i;
371 
372  globalResidual.data()[offset_u+stride_u*r_l2g.data()[eN_i]]+=elementResidual_u[i];
373  }//i
374  }//elements
375  //
376  //loop over exterior element boundaries to calculate levelset gradient
377  //
378  //ebNE is the Exterior element boundary INdex
379  //ebN is the element boundary INdex
380  //eN is the element index
381  for (int ebNE = 0; ebNE < nExteriorElementBoundaries_global; ebNE++)
382  {
383  register int ebN = exteriorElementBoundariesArray.data()[ebNE],
384  eN = elementBoundaryElementsArray.data()[ebN*2+0],
385  ebN_local = elementBoundaryLocalElementBoundariesArray.data()[ebN*2+0];
386  //eN_nDOF_trial_element = eN*nDOF_trial_element;
387  //register double elementResidual_u[nDOF_test_element];
388  double element_u[nDOF_trial_element];
389  for (int i=0;i<nDOF_test_element;i++)
390  {
391  register int eN_i=eN*nDOF_test_element+i;
392  element_u[i] = u_dof.data()[u_l2g.data()[eN_i]];
393  }//i
394  for (int kb=0;kb<nQuadraturePoints_elementBoundary;kb++)
395  {
396  register int ebNE_kb = ebNE*nQuadraturePoints_elementBoundary+kb,
397  ebNE_kb_nSpace = ebNE_kb*nSpace,
398  ebN_local_kb = ebN_local*nQuadraturePoints_elementBoundary+kb,
399  ebN_local_kb_nSpace = ebN_local_kb*nSpace;
400  register double u_ext=0.0,
401  grad_u_ext[nSpace],
402  //m_ext=0.0,
403  //dm_ext=0.0,
404  //H_ext=0.0,
405  //dH_ext[nSpace],
406  //flux_ext=0.0,
407  //bc_u_ext=0.0,
408  //bc_grad_u_ext[nSpace],
409  //bc_m_ext=0.0,
410  //bc_dm_ext=0.0,
411  //bc_H_ext=0.0,
412  //bc_dH_ext[nSpace],
413  jac_ext[nSpace*nSpace],
414  jacDet_ext,
415  jacInv_ext[nSpace*nSpace],
416  boundaryJac[nSpace*(nSpace-1)],
417  metricTensor[(nSpace-1)*(nSpace-1)],
418  metricTensorDetSqrt,
419  dS,
420  //u_test_dS[nDOF_test_element],
421  u_grad_trial_trace[nDOF_trial_element*nSpace],
422  normal[nSpace],x_ext,y_ext,z_ext,
423  G[nSpace*nSpace],G_dd_G,tr_G,norm;
424  //
425  //calculate the solution and gradients at quadrature points
426  //
427  ck.calculateMapping_elementBoundary(eN,
428  ebN_local,
429  kb,
430  ebN_local_kb,
431  mesh_dof.data(),
432  mesh_l2g.data(),
433  mesh_trial_trace_ref.data(),
434  mesh_grad_trial_trace_ref.data(),
435  boundaryJac_ref.data(),
436  jac_ext,
437  jacDet_ext,
438  jacInv_ext,
439  boundaryJac,
440  metricTensor,
441  metricTensorDetSqrt,
442  normal_ref.data(),
443  normal,
444  x_ext,y_ext,z_ext);
445  dS = metricTensorDetSqrt*dS_ref.data()[kb];
446  //get the metric tensor
447  //cek todo use symmetry
448  ck.calculateG(jacInv_ext,G,G_dd_G,tr_G);
449  //compute shape and solution information
450  //shape
451  ck.gradTrialFromRef(&u_grad_trial_trace_ref.data()[ebN_local_kb_nSpace*nDOF_trial_element],jacInv_ext,u_grad_trial_trace);
452  //solution and gradients
453  ck.valFromElementDOF(element_u,&u_trial_trace_ref.data()[ebN_local_kb*nDOF_test_element],u_ext);
454  ck.gradFromElementDOF(element_u,u_grad_trial_trace,grad_u_ext);
455 
456  ebqe_u.data()[ebNE_kb] = u_ext;
457  norm = 1.0e-8;
458  for (int I=0;I<nSpace;I++)
459  norm += grad_u_ext[I]*grad_u_ext[I];
460  norm = sqrt(norm);
461  for (int I=0;I<nSpace;I++)
462  ebqe_n.data()[ebNE_kb_nSpace+I] = grad_u_ext[I]/norm;
463  }//kb
464  }//ebNE
465  }
466 
467  inline void calculateElementJacobian(//element
468  double* mesh_trial_ref,
469  double* mesh_grad_trial_ref,
470  double* mesh_dof,
471  int* mesh_l2g,
472  double* dV_ref,
473  double* u_trial_ref,
474  double* u_grad_trial_ref,
475  double* u_test_ref,
476  double* u_grad_test_ref,
477  //element boundary
478  double* mesh_trial_trace_ref,
479  double* mesh_grad_trial_trace_ref,
480  double* dS_ref,
481  double* u_trial_trace_ref,
482  double* u_grad_trial_trace_ref,
483  double* u_test_trace_ref,
484  double* u_grad_test_trace_ref,
485  double* normal_ref,
486  double* boundaryJac_ref,
487  //physics
488  int nElements_global,
489  double useMetrics,
490  double epsFactHeaviside,
491  double epsFactDirac,
492  double epsFactDiffusion,
493  int* u_l2g,
494  double* elementDiameter,
495  double* nodeDiametersArray,
496  double* u_dof,
497  // double* u_trial,
498  // double* u_grad_trial,
499  // double* u_test_dV,
500  // double* u_grad_test_dV,
501  double* q_phi,
502  double* q_normal_phi,
503  double* q_H,
504  double* q_porosity,
505  double* elementJacobian_u_u,
506  double* element_u,
507  int eN)
508  {
509  for (int i=0;i<nDOF_test_element;i++)
510  for (int j=0;j<nDOF_trial_element;j++)
511  {
512  elementJacobian_u_u[i*nDOF_trial_element+j]=0.0;
513  }
514  double epsHeaviside,epsDirac,epsDiffusion;
515  for (int k=0;k<nQuadraturePoints_element;k++)
516  {
517  int eN_k = eN*nQuadraturePoints_element+k, //index to a scalar at a quadrature point
518  eN_k_nSpace = eN_k*nSpace;
519  //eN_nDOF_trial_element = eN*nDOF_trial_element; //index to a vector at a quadrature point
520  gf.set_quad(k);
521  //declare local storage
522  register double u=0.0,
523  grad_u[nSpace],
524  r=0.0,dr=0.0,
525  jac[nSpace*nSpace],
526  jacDet,
527  jacInv[nSpace*nSpace],
528  u_grad_trial[nDOF_trial_element*nSpace],
529  dV,
530  u_test_dV[nDOF_test_element],
531  u_grad_test_dV[nDOF_test_element*nSpace],
532  x,y,z,
533  G[nSpace*nSpace],G_dd_G,tr_G,h_phi;
534  //
535  //calculate solution and gradients at quadrature points
536  //
537  ck.calculateMapping_element(eN,
538  k,
539  mesh_dof,
540  mesh_l2g,
541  mesh_trial_ref,
542  mesh_grad_trial_ref,
543  jac,
544  jacDet,
545  jacInv,
546  x,y,z);
547  ck.calculateH_element(eN,
548  k,
549  nodeDiametersArray,
550  mesh_l2g,
551  mesh_trial_ref,
552  h_phi);
553  //get the physical integration weight
554  dV = fabs(jacDet)*dV_ref[k];
555  ck.calculateG(jacInv,G,G_dd_G,tr_G);
556 
557  /* double dir[nSpace]; */
558  /* double norm = 1.0e-8; */
559  /* for (int I=0;I<nSpace;I++) */
560  /* norm += q_normal_phi[eN_k_nSpace+I]*q_normal_phi[eN_k_nSpace+I]; */
561  /* norm = sqrt(norm); */
562  /* for (int I=0;I<nSpace;I++) */
563  /* dir[I] = q_normal_phi[eN_k_nSpace+I]/norm; */
564  /* ck.calculateGScale(G,dir,h_phi); */
565 
566 
567  //get the trial function gradients
568  ck.gradTrialFromRef(&u_grad_trial_ref[k*nDOF_trial_element*nSpace],jacInv,u_grad_trial);
569  //get the solution
570  ck.valFromElementDOF(element_u,&u_trial_ref[k*nDOF_trial_element],u);
571  //get the solution gradients
572  ck.gradFromElementDOF(element_u,u_grad_trial,grad_u);
573  //precalculate test function products with integration weights
574  for (int j=0;j<nDOF_trial_element;j++)
575  {
576  u_test_dV[j] = u_test_ref[k*nDOF_trial_element+j]*dV;
577  for (int I=0;I<nSpace;I++)
578  {
579  u_grad_test_dV[j*nSpace+I] = u_grad_trial[j*nSpace+I]*dV;//cek warning won't work for Petrov-Galerkin
580  }
581  }
582  //
583  //calculate pde coefficients and derivatives at quadrature points
584  //
585  epsHeaviside=epsFactHeaviside*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter[eN]);
586  epsDirac =epsFactDirac* (useMetrics*h_phi+(1.0-useMetrics)*elementDiameter[eN]);
587  epsDiffusion=epsFactDiffusion*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter[eN]);
588  // *(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter[eN]);
589  evaluateCoefficients(epsHeaviside,
590  epsDirac,
591  q_phi[eN_k],
592  q_H[eN_k],
593  u,
594  q_porosity[eN_k],
595  r,
596  dr);
597  for(int i=0;i<nDOF_test_element;i++)
598  {
599  //int eN_k_i=eN_k*nDOF_test_element+i;
600  //int eN_k_i_nSpace=eN_k_i*nSpace;
601  int i_nSpace=i*nSpace;
602  for(int j=0;j<nDOF_trial_element;j++)
603  {
604  //int eN_k_j=eN_k*nDOF_trial_element+j;
605  //int eN_k_j_nSpace = eN_k_j*nSpace;
606  int j_nSpace = j*nSpace;
607  elementJacobian_u_u[i*nDOF_trial_element+j] +=
608  ck.ReactionJacobian_weak(dr,u_trial_ref[k*nDOF_trial_element+j],u_test_dV[i]) +
609  ck.NumericalDiffusionJacobian(epsDiffusion,&u_grad_trial[j_nSpace],&u_grad_test_dV[i_nSpace]);
610  }//j
611  }//i
612  }//k
613  }
614 
615  void calculateJacobian(arguments_dict& args,
616  bool useExact)
617  {
618  xt::pyarray<double>& mesh_trial_ref = args.array<double>("mesh_trial_ref");
619  xt::pyarray<double>& mesh_grad_trial_ref = args.array<double>("mesh_grad_trial_ref");
620  xt::pyarray<double>& mesh_dof = args.array<double>("mesh_dof");
621  xt::pyarray<int>& mesh_l2g = args.array<int>("mesh_l2g");
622  xt::pyarray<double>& x_ref = args.array<double>("x_ref");
623  xt::pyarray<double>& dV_ref = args.array<double>("dV_ref");
624  xt::pyarray<double>& u_trial_ref = args.array<double>("u_trial_ref");
625  xt::pyarray<double>& u_grad_trial_ref = args.array<double>("u_grad_trial_ref");
626  xt::pyarray<double>& u_test_ref = args.array<double>("u_test_ref");
627  xt::pyarray<double>& u_grad_test_ref = args.array<double>("u_grad_test_ref");
628  xt::pyarray<double>& mesh_trial_trace_ref = args.array<double>("mesh_trial_trace_ref");
629  xt::pyarray<double>& mesh_grad_trial_trace_ref = args.array<double>("mesh_grad_trial_trace_ref");
630  xt::pyarray<double>& dS_ref = args.array<double>("dS_ref");
631  xt::pyarray<double>& u_trial_trace_ref = args.array<double>("u_trial_trace_ref");
632  xt::pyarray<double>& u_grad_trial_trace_ref = args.array<double>("u_grad_trial_trace_ref");
633  xt::pyarray<double>& u_test_trace_ref = args.array<double>("u_test_trace_ref");
634  xt::pyarray<double>& u_grad_test_trace_ref = args.array<double>("u_grad_test_trace_ref");
635  xt::pyarray<double>& normal_ref = args.array<double>("normal_ref");
636  xt::pyarray<double>& boundaryJac_ref = args.array<double>("boundaryJac_ref");
637  int nElements_global = args.scalar<int>("nElements_global");
638  double useMetrics = args.scalar<double>("useMetrics");
639  double epsFactHeaviside = args.scalar<double>("epsFactHeaviside");
640  double epsFactDirac = args.scalar<double>("epsFactDirac");
641  double epsFactDiffusion = args.scalar<double>("epsFactDiffusion");
642  xt::pyarray<int>& u_l2g = args.array<int>("u_l2g");
643  xt::pyarray<double>& elementDiameter = args.array<double>("elementDiameter");
644  xt::pyarray<double>& nodeDiametersArray = args.array<double>("nodeDiametersArray");
645  xt::pyarray<double>& u_dof = args.array<double>("u_dof");
646  xt::pyarray<double>& phi_dof = args.array<double>("phi_dof");
647  xt::pyarray<double>& q_phi = args.array<double>("q_phi");
648  xt::pyarray<double>& q_normal_phi = args.array<double>("q_normal_phi");
649  xt::pyarray<double>& q_H = args.array<double>("q_H");
650  xt::pyarray<double>& q_porosity = args.array<double>("q_porosity");
651  xt::pyarray<int>& csrRowIndeces_u_u = args.array<int>("csrRowIndeces_u_u");
652  xt::pyarray<int>& csrColumnOffsets_u_u = args.array<int>("csrColumnOffsets_u_u");
653  xt::pyarray<double>& globalJacobian = args.array<double>("globalJacobian");
654  //
655  //loop over elements to compute volume integrals and load them into the element Jacobians and global Jacobian
656  //
657  gf.useExact = useExact;
658  for(int eN=0;eN<nElements_global;eN++)
659  {
660  register double elementJacobian_u_u[nDOF_test_element*nDOF_trial_element],element_u[nDOF_trial_element],element_phi[nDOF_trial_element];
661  for (int j=0;j<nDOF_trial_element;j++)
662  {
663  register int eN_j = eN*nDOF_trial_element+j;
664  element_u[j] = u_dof.data()[u_l2g.data()[eN_j]];
665  element_phi[j] = phi_dof.data()[u_l2g.data()[eN_j]] + element_u[j];
666  }
667  double element_nodes[nDOF_mesh_trial_element*3];
668  for (int i=0;i<nDOF_mesh_trial_element;i++)
669  {
670  register int eN_i=eN*nDOF_mesh_trial_element+i;
671  for(int I=0;I<3;I++)
672  element_nodes[i*3 + I] = mesh_dof.data()[mesh_l2g.data()[eN_i]*3 + I];
673  }//i
674  gf.calculate(element_phi, element_nodes, x_ref.data(),false);
675  calculateElementJacobian(mesh_trial_ref.data(),
676  mesh_grad_trial_ref.data(),
677  mesh_dof.data(),
678  mesh_l2g.data(),
679  dV_ref.data(),
680  u_trial_ref.data(),
681  u_grad_trial_ref.data(),
682  u_test_ref.data(),
683  u_grad_test_ref.data(),
684  mesh_trial_trace_ref.data(),
685  mesh_grad_trial_trace_ref.data(),
686  dS_ref.data(),
687  u_trial_trace_ref.data(),
688  u_grad_trial_trace_ref.data(),
689  u_test_trace_ref.data(),
690  u_grad_test_trace_ref.data(),
691  normal_ref.data(),
692  boundaryJac_ref.data(),
693  nElements_global,
694  useMetrics,
695  epsFactHeaviside,
696  epsFactDirac,
697  epsFactDiffusion,
698  u_l2g.data(),
699  elementDiameter.data(),
700  nodeDiametersArray.data(),
701  u_dof.data(),
702  q_phi.data(),
703  q_normal_phi.data(),
704  q_H.data(),
705  q_porosity.data(),
706  elementJacobian_u_u,
707  element_u,
708  eN);
709  //
710  //load into element Jacobian into global Jacobian
711  //
712  for (int i=0;i<nDOF_test_element;i++)
713  {
714  int eN_i = eN*nDOF_test_element+i;
715  for (int j=0;j<nDOF_trial_element;j++)
716  {
717  int eN_i_j = eN_i*nDOF_trial_element+j;
718 
719  globalJacobian.data()[csrRowIndeces_u_u.data()[eN_i] + csrColumnOffsets_u_u.data()[eN_i_j]] += elementJacobian_u_u[i*nDOF_trial_element+j];
720  }//j
721  }//i
722  }//elements
723  }//computeJacobian
724  void elementSolve(arguments_dict& args)
725  {
726  xt::pyarray<double>& mesh_trial_ref = args.array<double>("mesh_trial_ref");
727  xt::pyarray<double>& mesh_grad_trial_ref = args.array<double>("mesh_grad_trial_ref");
728  xt::pyarray<double>& mesh_dof = args.array<double>("mesh_dof");
729  xt::pyarray<int>& mesh_l2g = args.array<int>("mesh_l2g");
730  xt::pyarray<double>& dV_ref = args.array<double>("dV_ref");
731  xt::pyarray<double>& u_trial_ref = args.array<double>("u_trial_ref");
732  xt::pyarray<double>& u_grad_trial_ref = args.array<double>("u_grad_trial_ref");
733  xt::pyarray<double>& u_test_ref = args.array<double>("u_test_ref");
734  xt::pyarray<double>& u_grad_test_ref = args.array<double>("u_grad_test_ref");
735  xt::pyarray<double>& mesh_trial_trace_ref = args.array<double>("mesh_trial_trace_ref");
736  xt::pyarray<double>& mesh_grad_trial_trace_ref = args.array<double>("mesh_grad_trial_trace_ref");
737  xt::pyarray<double>& dS_ref = args.array<double>("dS_ref");
738  xt::pyarray<double>& u_trial_trace_ref = args.array<double>("u_trial_trace_ref");
739  xt::pyarray<double>& u_grad_trial_trace_ref = args.array<double>("u_grad_trial_trace_ref");
740  xt::pyarray<double>& u_test_trace_ref = args.array<double>("u_test_trace_ref");
741  xt::pyarray<double>& u_grad_test_trace_ref = args.array<double>("u_grad_test_trace_ref");
742  xt::pyarray<double>& normal_ref = args.array<double>("normal_ref");
743  xt::pyarray<double>& boundaryJac_ref = args.array<double>("boundaryJac_ref");
744  int nElements_global = args.scalar<int>("nElements_global");
745  double useMetrics = args.scalar<double>("useMetrics");
746  double epsFactHeaviside = args.scalar<double>("epsFactHeaviside");
747  double epsFactDirac = args.scalar<double>("epsFactDirac");
748  double epsFactDiffusion = args.scalar<double>("epsFactDiffusion");
749  xt::pyarray<int>& u_l2g = args.array<int>("u_l2g");
750  xt::pyarray<double>& elementDiameter = args.array<double>("elementDiameter");
751  xt::pyarray<double>& nodeDiametersArray = args.array<double>("nodeDiametersArray");
752  xt::pyarray<double>& u_dof = args.array<double>("u_dof");
753  xt::pyarray<double>& q_phi = args.array<double>("q_phi");
754  xt::pyarray<double>& q_normal_phi = args.array<double>("q_normal_phi");
755  xt::pyarray<double>& ebqe_phi = args.array<double>("ebqe_phi");
756  xt::pyarray<double>& ebqe_normal_phi = args.array<double>("ebqe_normal_phi");
757  xt::pyarray<double>& q_H = args.array<double>("q_H");
758  xt::pyarray<double>& q_u = args.array<double>("q_u");
759  xt::pyarray<double>& q_n = args.array<double>("q_n");
760  xt::pyarray<double>& ebqe_u = args.array<double>("ebqe_u");
761  xt::pyarray<double>& ebqe_n = args.array<double>("ebqe_n");
762  xt::pyarray<double>& q_r = args.array<double>("q_r");
763  xt::pyarray<double>& q_porosity = args.array<double>("q_porosity");
764  int offset_u = args.scalar<int>("offset_u");
765  int stride_u = args.scalar<int>("stride_u");
766  xt::pyarray<double>& globalResidual = args.array<double>("globalResidual");
767  int nExteriorElementBoundaries_global = args.scalar<int>("nExteriorElementBoundaries_global");
768  xt::pyarray<int>& exteriorElementBoundariesArray = args.array<int>("exteriorElementBoundariesArray");
769  xt::pyarray<int>& elementBoundaryElementsArray = args.array<int>("elementBoundaryElementsArray");
770  xt::pyarray<int>& elementBoundaryLocalElementBoundariesArray = args.array<int>("elementBoundaryLocalElementBoundariesArray");
771  int maxIts = args.scalar<int>("maxIts");
772  double atol = args.scalar<double>("atol");
773  //
774  //loop over elements to compute volume integrals and load them into element and global residual
775  //
776  //eN is the element index
777  //eN_k is the quadrature point index for a scalar
778  //eN_k_nSpace is the quadrature point index for a vector
779  //eN_i is the element test function index
780  //eN_j is the element trial function index
781  //eN_k_j is the quadrature point index for a trial function
782  //eN_k_i is the quadrature point index for a trial function
783  for(int eN=0;eN<nElements_global;eN++)
784  {
785  //declare local storage for element residual and initialize
786  register double element_u[nDOF_test_element],
787  element_du[nDOF_test_element],
788  elementResidual_u[nDOF_test_element],
789  elementJacobian_u_u[nDOF_test_element*nDOF_trial_element],scale=1.0;
790  register PROTEUS_LAPACK_INTEGER elementPivots[nDOF_test_element],
791  elementColPivots[nDOF_test_element];
792  //double epsHeaviside,epsDirac,epsDiffusion;
793  for (int i=0;i<nDOF_test_element;i++)
794  {
795  element_u[i]=0.0;
796  }//i
797  calculateElementResidual(mesh_trial_ref.data(),
798  mesh_grad_trial_ref.data(),
799  mesh_dof.data(),
800  mesh_l2g.data(),
801  dV_ref.data(),
802  u_trial_ref.data(),
803  u_grad_trial_ref.data(),
804  u_test_ref.data(),
805  u_grad_test_ref.data(),
806  mesh_trial_trace_ref.data(),
807  mesh_grad_trial_trace_ref.data(),
808  dS_ref.data(),
809  u_trial_trace_ref.data(),
810  u_grad_trial_trace_ref.data(),
811  u_test_trace_ref.data(),
812  u_grad_test_trace_ref.data(),
813  normal_ref.data(),
814  boundaryJac_ref.data(),
815  nElements_global,
816  useMetrics,
817  epsFactHeaviside,
818  epsFactDirac,
819  epsFactDiffusion,
820  u_l2g.data(),
821  elementDiameter.data(),
822  nodeDiametersArray.data(),
823  u_dof.data(),
824  q_phi.data(),
825  q_normal_phi.data(),
826  ebqe_phi.data(),
827  ebqe_normal_phi.data(),
828  q_H.data(),
829  q_u.data(),
830  q_n.data(),
831  ebqe_u.data(),
832  ebqe_n.data(),
833  q_r.data(),
834  q_porosity.data(),
835  offset_u,stride_u,
836  elementResidual_u,
837  nExteriorElementBoundaries_global,
838  exteriorElementBoundariesArray.data(),
839  elementBoundaryElementsArray.data(),
840  elementBoundaryLocalElementBoundariesArray.data(),
841  element_u,
842  eN);
843  //compute l2 norm
844  double resNorm=0.0;
845  for (int i=0;i<nDOF_test_element;i++)
846  {
847  resNorm += elementResidual_u[i];
848  }//i
849  resNorm = fabs(resNorm);
850  //now do Newton
851  int its=0;
852  //std::cout<<"element "<<eN<<std::endl;
853  //std::cout<<"resNorm0 "<<resNorm<<std::endl;
854  while (resNorm >= atol && its < maxIts)
855  {
856  its+=1;
857  calculateElementJacobian(mesh_trial_ref.data(),
858  mesh_grad_trial_ref.data(),
859  mesh_dof.data(),
860  mesh_l2g.data(),
861  dV_ref.data(),
862  u_trial_ref.data(),
863  u_grad_trial_ref.data(),
864  u_test_ref.data(),
865  u_grad_test_ref.data(),
866  mesh_trial_trace_ref.data(),
867  mesh_grad_trial_trace_ref.data(),
868  dS_ref.data(),
869  u_trial_trace_ref.data(),
870  u_grad_trial_trace_ref.data(),
871  u_test_trace_ref.data(),
872  u_grad_test_trace_ref.data(),
873  normal_ref.data(),
874  boundaryJac_ref.data(),
875  nElements_global,
876  useMetrics,
877  epsFactHeaviside,
878  epsFactDirac,
879  epsFactDiffusion,
880  u_l2g.data(),
881  elementDiameter.data(),
882  nodeDiametersArray.data(),
883  u_dof.data(),
884  q_phi.data(),
885  q_normal_phi.data(),
886  q_H.data(),
887  q_porosity.data(),
888  elementJacobian_u_u,
889  element_u,
890  eN);
891  for (int i=0;i<nDOF_test_element;i++)
892  {
893  element_du[i] = -elementResidual_u[i];
894  elementPivots[i] = ((PROTEUS_LAPACK_INTEGER)0);
895  elementColPivots[i]=((PROTEUS_LAPACK_INTEGER)0);
896  /* std::cout<<"element jacobian"<<std::endl; */
897  /* for (int j=0;j<nDOF_test_element;j++) */
898  /* { */
899  /* std::cout<<elementJacobian_u_u[i*nDOF_trial_element+j]<<'\t'; */
900  /* } */
901  /* std::cout<<std::endl; */
902  }//i
903  //factor
904  PROTEUS_LAPACK_INTEGER La_N=((PROTEUS_LAPACK_INTEGER)nDOF_test_element),
905  INFO=0;
906  dgetc2_(&La_N,
907  elementJacobian_u_u,
908  &La_N,
909  elementPivots,
910  elementColPivots,
911  &INFO);
912  //solve
913  dgesc2_(&La_N,
914  elementJacobian_u_u,
915  &La_N,
916  element_du,
917  elementPivots,
918  elementColPivots,
919  &scale);
920  double resNormNew = resNorm,lambda=1.0;
921  int lsIts=0;
922  while (resNormNew > 0.99*resNorm && lsIts < 100)
923  {
924  //apply correction
925  for (int i=0;i<nDOF_test_element;i++)
926  {
927  element_u[i] += lambda*element_du[i];
928  }//i
929  lambda /= 2.0;
930  //compute new residual
931  calculateElementResidual(mesh_trial_ref.data(),
932  mesh_grad_trial_ref.data(),
933  mesh_dof.data(),
934  mesh_l2g.data(),
935  dV_ref.data(),
936  u_trial_ref.data(),
937  u_grad_trial_ref.data(),
938  u_test_ref.data(),
939  u_grad_test_ref.data(),
940  mesh_trial_trace_ref.data(),
941  mesh_grad_trial_trace_ref.data(),
942  dS_ref.data(),
943  u_trial_trace_ref.data(),
944  u_grad_trial_trace_ref.data(),
945  u_test_trace_ref.data(),
946  u_grad_test_trace_ref.data(),
947  normal_ref.data(),
948  boundaryJac_ref.data(),
949  nElements_global,
950  useMetrics,
951  epsFactHeaviside,
952  epsFactDirac,
953  epsFactDiffusion,
954  u_l2g.data(),
955  elementDiameter.data(),
956  nodeDiametersArray.data(),
957  u_dof.data(),
958  q_phi.data(),
959  q_normal_phi.data(),
960  ebqe_phi.data(),
961  ebqe_normal_phi.data(),
962  q_H.data(),
963  q_u.data(),
964  q_n.data(),
965  ebqe_u.data(),
966  ebqe_n.data(),
967  q_r.data(),
968  q_porosity.data(),
969  offset_u,stride_u,
970  elementResidual_u,
971  nExteriorElementBoundaries_global,
972  exteriorElementBoundariesArray.data(),
973  elementBoundaryElementsArray.data(),
974  elementBoundaryLocalElementBoundariesArray.data(),
975  element_u,
976  eN);
977  lsIts +=1;
978  //compute l2 norm
979  resNormNew=0.0;
980  for (int i=0;i<nDOF_test_element;i++)
981  {
982  resNormNew += elementResidual_u[i];
983  std::cout<<"element_u["<<i<<"] "<<element_u[i]<<std::endl;
984  std::cout<<"elementResidual_u["<<i<<"] "<<elementResidual_u[i]<<std::endl;
985  }//i
986  resNormNew = fabs(resNormNew);
987  }
988  resNorm = resNormNew;
989  std::cout<<"INFO "<<INFO<<std::endl;
990  std::cout<<"resNorm["<<its<<"] "<<resNorm<<std::endl;
991  }
992  }//elements
993  }
994  void elementConstantSolve(arguments_dict& args)
995  {
996  xt::pyarray<double>& mesh_trial_ref = args.array<double>("mesh_trial_ref");
997  xt::pyarray<double>& mesh_grad_trial_ref = args.array<double>("mesh_grad_trial_ref");
998  xt::pyarray<double>& mesh_dof = args.array<double>("mesh_dof");
999  xt::pyarray<int>& mesh_l2g = args.array<int>("mesh_l2g");
1000  xt::pyarray<double>& dV_ref = args.array<double>("dV_ref");
1001  xt::pyarray<double>& u_trial_ref = args.array<double>("u_trial_ref");
1002  xt::pyarray<double>& u_grad_trial_ref = args.array<double>("u_grad_trial_ref");
1003  xt::pyarray<double>& u_test_ref = args.array<double>("u_test_ref");
1004  xt::pyarray<double>& u_grad_test_ref = args.array<double>("u_grad_test_ref");
1005  xt::pyarray<double>& mesh_trial_trace_ref = args.array<double>("mesh_trial_trace_ref");
1006  xt::pyarray<double>& mesh_grad_trial_trace_ref = args.array<double>("mesh_grad_trial_trace_ref");
1007  xt::pyarray<double>& dS_ref = args.array<double>("dS_ref");
1008  xt::pyarray<double>& u_trial_trace_ref = args.array<double>("u_trial_trace_ref");
1009  xt::pyarray<double>& u_grad_trial_trace_ref = args.array<double>("u_grad_trial_trace_ref");
1010  xt::pyarray<double>& u_test_trace_ref = args.array<double>("u_test_trace_ref");
1011  xt::pyarray<double>& u_grad_test_trace_ref = args.array<double>("u_grad_test_trace_ref");
1012  xt::pyarray<double>& normal_ref = args.array<double>("normal_ref");
1013  xt::pyarray<double>& boundaryJac_ref = args.array<double>("boundaryJac_ref");
1014  int nElements_global = args.scalar<int>("nElements_global");
1015  double useMetrics = args.scalar<double>("useMetrics");
1016  double epsFactHeaviside = args.scalar<double>("epsFactHeaviside");
1017  double epsFactDirac = args.scalar<double>("epsFactDirac");
1018  double epsFactDiffusion = args.scalar<double>("epsFactDiffusion");
1019  xt::pyarray<int>& u_l2g = args.array<int>("u_l2g");
1020  xt::pyarray<double>& elementDiameter = args.array<double>("elementDiameter");
1021  xt::pyarray<double>& nodeDiametersArray = args.array<double>("nodeDiametersArray");
1022  xt::pyarray<double>& u_dof = args.array<double>("u_dof");
1023  xt::pyarray<double>& q_phi = args.array<double>("q_phi");
1024  xt::pyarray<double>& q_normal_phi = args.array<double>("q_normal_phi");
1025  xt::pyarray<double>& ebqe_phi = args.array<double>("ebqe_phi");
1026  xt::pyarray<double>& ebqe_normal_phi = args.array<double>("ebqe_normal_phi");
1027  xt::pyarray<double>& q_H = args.array<double>("q_H");
1028  xt::pyarray<double>& q_u = args.array<double>("q_u");
1029  xt::pyarray<double>& q_n = args.array<double>("q_n");
1030  xt::pyarray<double>& ebqe_u = args.array<double>("ebqe_u");
1031  xt::pyarray<double>& ebqe_n = args.array<double>("ebqe_n");
1032  xt::pyarray<double>& q_r = args.array<double>("q_r");
1033  xt::pyarray<double>& q_porosity = args.array<double>("q_porosity");
1034  int offset_u = args.scalar<int>("offset_u");
1035  int stride_u = args.scalar<int>("stride_u");
1036  xt::pyarray<double>& globalResidual = args.array<double>("globalResidual");
1037  int nExteriorElementBoundaries_global = args.scalar<int>("nExteriorElementBoundaries_global");
1038  xt::pyarray<int>& exteriorElementBoundariesArray = args.array<int>("exteriorElementBoundariesArray");
1039  xt::pyarray<int>& elementBoundaryElementsArray = args.array<int>("elementBoundaryElementsArray");
1040  xt::pyarray<int>& elementBoundaryLocalElementBoundariesArray = args.array<int>("elementBoundaryLocalElementBoundariesArray");
1041  int maxIts = args.scalar<int>("maxIts");
1042  double atol = args.scalar<double>("atol");
1043  for(int eN=0;eN<nElements_global;eN++)
1044  {
1045  //declare local storage for element residual and initialize
1046  register double element_u[nDOF_test_element],elementConstant_u,
1047  elementResidual_u[nDOF_test_element],elementConstantResidual,
1048  elementJacobian_u_u[nDOF_test_element*nDOF_trial_element],elementConstantJacobian,resNorm;
1049  elementConstant_u=0.0;
1050  for (int i=0;i<nDOF_test_element;i++)
1051  {
1052  element_u[i]=elementConstant_u;
1053  }//i
1054  calculateElementResidual(mesh_trial_ref.data(),
1055  mesh_grad_trial_ref.data(),
1056  mesh_dof.data(),
1057  mesh_l2g.data(),
1058  dV_ref.data(),
1059  u_trial_ref.data(),
1060  u_grad_trial_ref.data(),
1061  u_test_ref.data(),
1062  u_grad_test_ref.data(),
1063  mesh_trial_trace_ref.data(),
1064  mesh_grad_trial_trace_ref.data(),
1065  dS_ref.data(),
1066  u_trial_trace_ref.data(),
1067  u_grad_trial_trace_ref.data(),
1068  u_test_trace_ref.data(),
1069  u_grad_test_trace_ref.data(),
1070  normal_ref.data(),
1071  boundaryJac_ref.data(),
1072  nElements_global,
1073  useMetrics,
1074  epsFactHeaviside,
1075  epsFactDirac,
1076  epsFactDiffusion,
1077  u_l2g.data(),
1078  elementDiameter.data(),
1079  nodeDiametersArray.data(),
1080  u_dof.data(),
1081  q_phi.data(),
1082  q_normal_phi.data(),
1083  ebqe_phi.data(),
1084  ebqe_normal_phi.data(),
1085  q_H.data(),
1086  q_u.data(),
1087  q_n.data(),
1088  ebqe_u.data(),
1089  ebqe_n.data(),
1090  q_r.data(),
1091  q_porosity.data(),
1092  offset_u,stride_u,
1093  elementResidual_u,
1094  nExteriorElementBoundaries_global,
1095  exteriorElementBoundariesArray.data(),
1096  elementBoundaryElementsArray.data(),
1097  elementBoundaryLocalElementBoundariesArray.data(),
1098  element_u,
1099  eN);
1100  //compute l2 norm
1101  elementConstantResidual=0.0;
1102  for (int i=0;i<nDOF_test_element;i++)
1103  {
1104  elementConstantResidual += elementResidual_u[i];
1105  }//i
1106  resNorm = fabs(elementConstantResidual);
1107  //now do Newton
1108  int its=0;
1109  //std::cout<<"element "<<eN<<std::endl;
1110  //std::cout<<"resNorm0 "<<resNorm<<std::endl;
1111  while (resNorm >= atol && its < maxIts)
1112  {
1113  its+=1;
1114  calculateElementJacobian(mesh_trial_ref.data(),
1115  mesh_grad_trial_ref.data(),
1116  mesh_dof.data(),
1117  mesh_l2g.data(),
1118  dV_ref.data(),
1119  u_trial_ref.data(),
1120  u_grad_trial_ref.data(),
1121  u_test_ref.data(),
1122  u_grad_test_ref.data(),
1123  mesh_trial_trace_ref.data(),
1124  mesh_grad_trial_trace_ref.data(),
1125  dS_ref.data(),
1126  u_trial_trace_ref.data(),
1127  u_grad_trial_trace_ref.data(),
1128  u_test_trace_ref.data(),
1129  u_grad_test_trace_ref.data(),
1130  normal_ref.data(),
1131  boundaryJac_ref.data(),
1132  nElements_global,
1133  useMetrics,
1134  epsFactHeaviside,
1135  epsFactDirac,
1136  epsFactDiffusion,
1137  u_l2g.data(),
1138  elementDiameter.data(),
1139  nodeDiametersArray.data(),
1140  u_dof.data(),
1141  q_phi.data(),
1142  q_normal_phi.data(),
1143  q_H.data(),
1144  q_porosity.data(),
1145  elementJacobian_u_u,
1146  element_u,
1147  eN);
1148  elementConstantJacobian=0.0;
1149  for (int i=0;i<nDOF_test_element;i++)
1150  {
1151  for (int j=0;j<nDOF_test_element;j++)
1152  {
1153  elementConstantJacobian += elementJacobian_u_u[i*nDOF_trial_element+j];
1154  }
1155  }//i
1156  std::cout<<"elementConstantJacobian "<<elementConstantJacobian<<std::endl;
1157  //apply correction
1158  elementConstant_u -= elementConstantResidual/(elementConstantJacobian+1.0e-8);
1159  for (int i=0;i<nDOF_test_element;i++)
1160  {
1161  element_u[i] = elementConstant_u;
1162  }//i
1163  //compute new residual
1164  calculateElementResidual(mesh_trial_ref.data(),
1165  mesh_grad_trial_ref.data(),
1166  mesh_dof.data(),
1167  mesh_l2g.data(),
1168  dV_ref.data(),
1169  u_trial_ref.data(),
1170  u_grad_trial_ref.data(),
1171  u_test_ref.data(),
1172  u_grad_test_ref.data(),
1173  mesh_trial_trace_ref.data(),
1174  mesh_grad_trial_trace_ref.data(),
1175  dS_ref.data(),
1176  u_trial_trace_ref.data(),
1177  u_grad_trial_trace_ref.data(),
1178  u_test_trace_ref.data(),
1179  u_grad_test_trace_ref.data(),
1180  normal_ref.data(),
1181  boundaryJac_ref.data(),
1182  nElements_global,
1183  useMetrics,
1184  epsFactHeaviside,
1185  epsFactDirac,
1186  epsFactDiffusion,
1187  u_l2g.data(),
1188  elementDiameter.data(),
1189  nodeDiametersArray.data(),
1190  u_dof.data(),
1191  q_phi.data(),
1192  q_normal_phi.data(),
1193  ebqe_phi.data(),
1194  ebqe_normal_phi.data(),
1195  q_H.data(),
1196  q_u.data(),
1197  q_n.data(),
1198  ebqe_u.data(),
1199  ebqe_n.data(),
1200  q_r.data(),
1201  q_porosity.data(),
1202  offset_u,stride_u,
1203  elementResidual_u,
1204  nExteriorElementBoundaries_global,
1205  exteriorElementBoundariesArray.data(),
1206  elementBoundaryElementsArray.data(),
1207  elementBoundaryLocalElementBoundariesArray.data(),
1208  element_u,
1209  eN);
1210  //compute l2 norm
1211  elementConstantResidual=0.0;
1212  for (int i=0;i<nDOF_test_element;i++)
1213  {
1214  elementConstantResidual += elementResidual_u[i];
1215  }//i
1216  resNorm = fabs(elementConstantResidual);
1217  std::cout<<"resNorm["<<its<<"] "<<resNorm<<std::endl;
1218  }
1219  }//elements
1220  }
1221 
1222  std::tuple<double, double> globalConstantRJ(arguments_dict& args)
1223  {
1224  xt::pyarray<double>& mesh_trial_ref = args.array<double>("mesh_trial_ref");
1225  xt::pyarray<double>& mesh_grad_trial_ref = args.array<double>("mesh_grad_trial_ref");
1226  xt::pyarray<double>& mesh_dof = args.array<double>("mesh_dof");
1227  xt::pyarray<int>& mesh_l2g = args.array<int>("mesh_l2g");
1228  xt::pyarray<double>& dV_ref = args.array<double>("dV_ref");
1229  xt::pyarray<double>& u_trial_ref = args.array<double>("u_trial_ref");
1230  xt::pyarray<double>& u_grad_trial_ref = args.array<double>("u_grad_trial_ref");
1231  xt::pyarray<double>& u_test_ref = args.array<double>("u_test_ref");
1232  xt::pyarray<double>& u_grad_test_ref = args.array<double>("u_grad_test_ref");
1233  xt::pyarray<double>& mesh_trial_trace_ref = args.array<double>("mesh_trial_trace_ref");
1234  xt::pyarray<double>& mesh_grad_trial_trace_ref = args.array<double>("mesh_grad_trial_trace_ref");
1235  xt::pyarray<double>& dS_ref = args.array<double>("dS_ref");
1236  xt::pyarray<double>& u_trial_trace_ref = args.array<double>("u_trial_trace_ref");
1237  xt::pyarray<double>& u_grad_trial_trace_ref = args.array<double>("u_grad_trial_trace_ref");
1238  xt::pyarray<double>& u_test_trace_ref = args.array<double>("u_test_trace_ref");
1239  xt::pyarray<double>& u_grad_test_trace_ref = args.array<double>("u_grad_test_trace_ref");
1240  xt::pyarray<double>& normal_ref = args.array<double>("normal_ref");
1241  xt::pyarray<double>& boundaryJac_ref = args.array<double>("boundaryJac_ref");
1242  int nElements_owned = args.scalar<int>("nElements_owned");
1243  double useMetrics = args.scalar<double>("useMetrics");
1244  double epsFactHeaviside = args.scalar<double>("epsFactHeaviside");
1245  double epsFactDirac = args.scalar<double>("epsFactDirac");
1246  double epsFactDiffusion = args.scalar<double>("epsFactDiffusion");
1247  xt::pyarray<int>& u_l2g = args.array<int>("u_l2g");
1248  xt::pyarray<double>& elementDiameter = args.array<double>("elementDiameter");
1249  xt::pyarray<double>& nodeDiametersArray = args.array<double>("nodeDiametersArray");
1250  xt::pyarray<double>& u_dof = args.array<double>("u_dof");
1251  xt::pyarray<double>& q_phi = args.array<double>("q_phi");
1252  xt::pyarray<double>& q_normal_phi = args.array<double>("q_normal_phi");
1253  xt::pyarray<double>& ebqe_phi = args.array<double>("ebqe_phi");
1254  xt::pyarray<double>& ebqe_normal_phi = args.array<double>("ebqe_normal_phi");
1255  xt::pyarray<double>& q_H = args.array<double>("q_H");
1256  xt::pyarray<double>& q_u = args.array<double>("q_u");
1257  xt::pyarray<double>& q_n = args.array<double>("q_n");
1258  xt::pyarray<double>& ebqe_u = args.array<double>("ebqe_u");
1259  xt::pyarray<double>& ebqe_n = args.array<double>("ebqe_n");
1260  xt::pyarray<double>& q_r = args.array<double>("q_r");
1261  xt::pyarray<double>& q_porosity = args.array<double>("q_porosity");
1262  int offset_u = args.scalar<int>("offset_u");
1263  int stride_u = args.scalar<int>("stride_u");
1264  xt::pyarray<double>& globalResidual = args.array<double>("globalResidual");
1265  int nExteriorElementBoundaries_global = args.scalar<int>("nExteriorElementBoundaries_global");
1266  xt::pyarray<int>& exteriorElementBoundariesArray = args.array<int>("exteriorElementBoundariesArray");
1267  xt::pyarray<int>& elementBoundaryElementsArray = args.array<int>("elementBoundaryElementsArray");
1268  xt::pyarray<int>& elementBoundaryLocalElementBoundariesArray = args.array<int>("elementBoundaryLocalElementBoundariesArray");
1269  int maxIts = args.scalar<int>("maxIts");
1270  double atol = args.scalar<double>("atol");
1271  double constant_u = args.scalar<double>("constant_u");
1272  register double element_u[nDOF_test_element],
1273  elementResidual_u[nDOF_test_element],
1274  elementJacobian_u_u[nDOF_test_element*nDOF_trial_element];
1275  double constantResidual = 0.0;
1276  double constantJacobian = 0.0;
1277  for (int i=0;i<nDOF_trial_element;i++)
1278  {
1279  element_u[i]=constant_u;
1280  }//i
1281  //compute residual and Jacobian
1282  for(int eN=0;eN<nElements_owned;eN++)
1283  {
1284  calculateElementResidual(mesh_trial_ref.data(),
1285  mesh_grad_trial_ref.data(),
1286  mesh_dof.data(),
1287  mesh_l2g.data(),
1288  dV_ref.data(),
1289  u_trial_ref.data(),
1290  u_grad_trial_ref.data(),
1291  u_test_ref.data(),
1292  u_grad_test_ref.data(),
1293  mesh_trial_trace_ref.data(),
1294  mesh_grad_trial_trace_ref.data(),
1295  dS_ref.data(),
1296  u_trial_trace_ref.data(),
1297  u_grad_trial_trace_ref.data(),
1298  u_test_trace_ref.data(),
1299  u_grad_test_trace_ref.data(),
1300  normal_ref.data(),
1301  boundaryJac_ref.data(),
1302  nElements_owned,
1303  useMetrics,
1304  epsFactHeaviside,
1305  epsFactDirac,
1306  epsFactDiffusion,
1307  u_l2g.data(),
1308  elementDiameter.data(),
1309  nodeDiametersArray.data(),
1310  u_dof.data(),
1311  q_phi.data(),
1312  q_normal_phi.data(),
1313  ebqe_phi.data(),
1314  ebqe_normal_phi.data(),
1315  q_H.data(),
1316  q_u.data(),
1317  q_n.data(),
1318  ebqe_u.data(),
1319  ebqe_n.data(),
1320  q_r.data(),
1321  q_porosity.data(),
1322  offset_u,stride_u,
1323  elementResidual_u,
1324  nExteriorElementBoundaries_global,
1325  exteriorElementBoundariesArray.data(),
1326  elementBoundaryElementsArray.data(),
1327  elementBoundaryLocalElementBoundariesArray.data(),
1328  element_u,
1329  eN);
1330  //compute l2 norm
1331  for (int i=0;i<nDOF_test_element;i++)
1332  {
1333  constantResidual += elementResidual_u[i];
1334  }//i
1335  calculateElementJacobian(mesh_trial_ref.data(),
1336  mesh_grad_trial_ref.data(),
1337  mesh_dof.data(),
1338  mesh_l2g.data(),
1339  dV_ref.data(),
1340  u_trial_ref.data(),
1341  u_grad_trial_ref.data(),
1342  u_test_ref.data(),
1343  u_grad_test_ref.data(),
1344  mesh_trial_trace_ref.data(),
1345  mesh_grad_trial_trace_ref.data(),
1346  dS_ref.data(),
1347  u_trial_trace_ref.data(),
1348  u_grad_trial_trace_ref.data(),
1349  u_test_trace_ref.data(),
1350  u_grad_test_trace_ref.data(),
1351  normal_ref.data(),
1352  boundaryJac_ref.data(),
1353  nElements_owned,
1354  useMetrics,
1355  epsFactHeaviside,
1356  epsFactDirac,
1357  epsFactDiffusion,
1358  u_l2g.data(),
1359  elementDiameter.data(),
1360  nodeDiametersArray.data(),
1361  u_dof.data(),
1362  q_phi.data(),
1363  q_normal_phi.data(),
1364  q_H.data(),
1365  q_porosity.data(),
1366  elementJacobian_u_u,
1367  element_u,
1368  eN);
1369  for (int i=0;i<nDOF_test_element;i++)
1370  {
1371  for (int j=0;j<nDOF_test_element;j++)
1372  {
1373  constantJacobian += elementJacobian_u_u[i*nDOF_trial_element+j];
1374  }
1375  }//i
1376  }
1377  return std::tuple<double, double>(constantResidual, constantJacobian);
1378  }
1379 
1380  double calculateMass(arguments_dict& args,
1381  bool useExact)
1382  {
1383  xt::pyarray<double>& mesh_trial_ref = args.array<double>("mesh_trial_ref");
1384  xt::pyarray<double>& mesh_grad_trial_ref = args.array<double>("mesh_grad_trial_ref");
1385  xt::pyarray<double>& mesh_dof = args.array<double>("mesh_dof");
1386  xt::pyarray<int>& mesh_l2g = args.array<int>("mesh_l2g");
1387  xt::pyarray<double>& x_ref = args.array<double>("x_ref");
1388  xt::pyarray<double>& dV_ref = args.array<double>("dV_ref");
1389  xt::pyarray<double>& u_trial_ref = args.array<double>("u_trial_ref");
1390  xt::pyarray<double>& u_grad_trial_ref = args.array<double>("u_grad_trial_ref");
1391  xt::pyarray<double>& u_test_ref = args.array<double>("u_test_ref");
1392  xt::pyarray<double>& u_grad_test_ref = args.array<double>("u_grad_test_ref");
1393  xt::pyarray<double>& mesh_trial_trace_ref = args.array<double>("mesh_trial_trace_ref");
1394  xt::pyarray<double>& mesh_grad_trial_trace_ref = args.array<double>("mesh_grad_trial_trace_ref");
1395  xt::pyarray<double>& dS_ref = args.array<double>("dS_ref");
1396  xt::pyarray<double>& u_trial_trace_ref = args.array<double>("u_trial_trace_ref");
1397  xt::pyarray<double>& u_grad_trial_trace_ref = args.array<double>("u_grad_trial_trace_ref");
1398  xt::pyarray<double>& u_test_trace_ref = args.array<double>("u_test_trace_ref");
1399  xt::pyarray<double>& u_grad_test_trace_ref = args.array<double>("u_grad_test_trace_ref");
1400  xt::pyarray<double>& normal_ref = args.array<double>("normal_ref");
1401  xt::pyarray<double>& boundaryJac_ref = args.array<double>("boundaryJac_ref");
1402  int nElements_owned = args.scalar<int>("nElements_owned");
1403  double useMetrics = args.scalar<double>("useMetrics");
1404  double epsFactHeaviside = args.scalar<double>("epsFactHeaviside");
1405  double epsFactDirac = args.scalar<double>("epsFactDirac");
1406  double epsFactDiffusion = args.scalar<double>("epsFactDiffusion");
1407  xt::pyarray<int>& u_l2g = args.array<int>("u_l2g");
1408  xt::pyarray<double>& elementDiameter = args.array<double>("elementDiameter");
1409  xt::pyarray<double>& nodeDiametersArray = args.array<double>("nodeDiametersArray");
1410  xt::pyarray<double>& u_dof = args.array<double>("u_dof");
1411  xt::pyarray<double>& phi_dof = args.array<double>("phi_dof");
1412  xt::pyarray<double>& q_phi = args.array<double>("q_phi");
1413  xt::pyarray<double>& q_normal_phi = args.array<double>("q_normal_phi");
1414  xt::pyarray<double>& ebqe_phi = args.array<double>("ebqe_phi");
1415  xt::pyarray<double>& ebqe_normal_phi = args.array<double>("ebqe_normal_phi");
1416  xt::pyarray<double>& q_H = args.array<double>("q_H");
1417  xt::pyarray<double>& q_u = args.array<double>("q_u");
1418  xt::pyarray<double>& q_n = args.array<double>("q_n");
1419  xt::pyarray<double>& ebqe_u = args.array<double>("ebqe_u");
1420  xt::pyarray<double>& ebqe_n = args.array<double>("ebqe_n");
1421  xt::pyarray<double>& q_r = args.array<double>("q_r");
1422  xt::pyarray<double>& q_porosity = args.array<double>("q_porosity");
1423  int offset_u = args.scalar<int>("offset_u");
1424  int stride_u = args.scalar<int>("stride_u");
1425  xt::pyarray<double>& globalResidual = args.array<double>("globalResidual");
1426  int nExteriorElementBoundaries_global = args.scalar<int>("nExteriorElementBoundaries_global");
1427  xt::pyarray<int>& exteriorElementBoundariesArray = args.array<int>("exteriorElementBoundariesArray");
1428  xt::pyarray<int>& elementBoundaryElementsArray = args.array<int>("elementBoundaryElementsArray");
1429  xt::pyarray<int>& elementBoundaryLocalElementBoundariesArray = args.array<int>("elementBoundaryLocalElementBoundariesArray");
1430  double globalMass = 0.0;
1431  gf.useExact=useExact;
1432  for(int eN=0;eN<nElements_owned;eN++)
1433  {
1434  double epsHeaviside;
1435  //loop over quadrature points and compute integrands
1436  //declare local storage for element residual and initialize
1437  register double element_phi[nDOF_trial_element];
1438  for (int i=0;i<nDOF_test_element;i++)
1439  {
1440  register int eN_i=eN*nDOF_test_element+i;
1441  element_phi[i] = phi_dof.data()[u_l2g.data()[eN_i]];
1442  }//i
1443  double element_nodes[nDOF_mesh_trial_element*3];
1444  for (int i=0;i<nDOF_mesh_trial_element;i++)
1445  {
1446  register int eN_i=eN*nDOF_mesh_trial_element+i;
1447  for(int I=0;I<3;I++)
1448  element_nodes[i*3 + I] = mesh_dof.data()[mesh_l2g.data()[eN_i]*3 + I];
1449  }//i
1450  gf.calculate(element_phi, element_nodes, x_ref.data(),false);
1451  for (int k=0;k<nQuadraturePoints_element;k++)
1452  {
1453  //compute indeces and declare local storage
1454  register int eN_k = eN*nQuadraturePoints_element+k,
1455  eN_k_nSpace = eN_k*nSpace;
1456  //eN_nDOF_trial_element = eN*nDOF_trial_element;
1457  //double u=0.0,grad_u[nSpace],r=0.0,dr=0.0;
1458  double jac[nSpace*nSpace],
1459  jacDet,
1460  jacInv[nSpace*nSpace],
1461  //u_grad_trial[nDOF_trial_element*nSpace],
1462  //u_test_dV[nDOF_trial_element],
1463  //u_grad_test_dV[nDOF_test_element*nSpace],
1464  dV,x,y,z,
1465  G[nSpace*nSpace],G_dd_G,tr_G,h_phi;
1466  gf.set_quad(k);
1467  //
1468  //compute solution and gradients at quadrature points
1469  //
1470  ck.calculateMapping_element(eN,
1471  k,
1472  mesh_dof.data(),
1473  mesh_l2g.data(),
1474  mesh_trial_ref.data(),
1475  mesh_grad_trial_ref.data(),
1476  jac,
1477  jacDet,
1478  jacInv,
1479  x,y,z);
1480  ck.calculateH_element(eN,
1481  k,
1482  nodeDiametersArray.data(),
1483  mesh_l2g.data(),
1484  mesh_trial_ref.data(),
1485  h_phi);
1486  //get the physical integration weight
1487  dV = fabs(jacDet)*dV_ref.data()[k];
1488  ck.calculateG(jacInv,G,G_dd_G,tr_G);
1489  /* double dir[nSpace]; */
1490  /* double norm = 1.0e-8; */
1491  /* for (int I=0;I<nSpace;I++) */
1492  /* norm += q_normal_phi.data()[eN_k_nSpace+I]*q_normal_phi.data()[eN_k_nSpace+I]; */
1493  /* norm = sqrt(norm); */
1494  /* for (int I=0;I<nSpace;I++) */
1495  /* dir[I] = q_normal_phi.data()[eN_k_nSpace+I]/norm; */
1496  /* ck.calculateGScale(G,dir,h_phi); */
1497  epsHeaviside=epsFactHeaviside*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter.data()[eN]);
1498  globalMass += q_porosity[eN_k]*gf.H(epsHeaviside,q_phi.data()[eN_k])*dV;
1499  }//k
1500  }//elements
1501  return globalMass;
1502  }
1503 
1504  void setMassQuadrature(arguments_dict& args,
1505  bool useExact)
1506  {
1507  xt::pyarray<double>& mesh_trial_ref = args.array<double>("mesh_trial_ref");
1508  xt::pyarray<double>& mesh_grad_trial_ref = args.array<double>("mesh_grad_trial_ref");
1509  xt::pyarray<double>& mesh_dof = args.array<double>("mesh_dof");
1510  xt::pyarray<int>& mesh_l2g = args.array<int>("mesh_l2g");
1511  xt::pyarray<double>& x_ref = args.array<double>("x_ref");
1512  xt::pyarray<double>& dV_ref = args.array<double>("dV_ref");
1513  xt::pyarray<double>& u_trial_ref = args.array<double>("u_trial_ref");
1514  xt::pyarray<double>& u_grad_trial_ref = args.array<double>("u_grad_trial_ref");
1515  xt::pyarray<double>& u_test_ref = args.array<double>("u_test_ref");
1516  xt::pyarray<double>& u_grad_test_ref = args.array<double>("u_grad_test_ref");
1517  xt::pyarray<double>& mesh_trial_trace_ref = args.array<double>("mesh_trial_trace_ref");
1518  xt::pyarray<double>& mesh_grad_trial_trace_ref = args.array<double>("mesh_grad_trial_trace_ref");
1519  xt::pyarray<double>& dS_ref = args.array<double>("dS_ref");
1520  xt::pyarray<double>& u_trial_trace_ref = args.array<double>("u_trial_trace_ref");
1521  xt::pyarray<double>& u_grad_trial_trace_ref = args.array<double>("u_grad_trial_trace_ref");
1522  xt::pyarray<double>& u_test_trace_ref = args.array<double>("u_test_trace_ref");
1523  xt::pyarray<double>& u_grad_test_trace_ref = args.array<double>("u_grad_test_trace_ref");
1524  xt::pyarray<double>& normal_ref = args.array<double>("normal_ref");
1525  xt::pyarray<double>& boundaryJac_ref = args.array<double>("boundaryJac_ref");
1526  int nElements_global = args.scalar<int>("nElements_global");
1527  double useMetrics = args.scalar<double>("useMetrics");
1528  double epsFactHeaviside = args.scalar<double>("epsFactHeaviside");
1529  double epsFactDirac = args.scalar<double>("epsFactDirac");
1530  double epsFactDiffusion = args.scalar<double>("epsFactDiffusion");
1531  xt::pyarray<int>& phi_l2g = args.array<int>("phi_l2g");
1532  xt::pyarray<double>& elementDiameter = args.array<double>("elementDiameter");
1533  xt::pyarray<double>& nodeDiametersArray = args.array<double>("nodeDiametersArray");
1534  xt::pyarray<double>& phi_dof = args.array<double>("phi_dof");
1535  xt::pyarray<double>& q_phi = args.array<double>("q_phi");
1536  xt::pyarray<double>& q_normal_phi = args.array<double>("q_normal_phi");
1537  xt::pyarray<double>& ebqe_phi = args.array<double>("ebqe_phi");
1538  xt::pyarray<double>& ebqe_normal_phi = args.array<double>("ebqe_normal_phi");
1539  xt::pyarray<double>& q_H = args.array<double>("q_H");
1540  xt::pyarray<double>& q_u = args.array<double>("q_u");
1541  xt::pyarray<double>& q_n = args.array<double>("q_n");
1542  xt::pyarray<double>& ebqe_u = args.array<double>("ebqe_u");
1543  xt::pyarray<double>& ebqe_n = args.array<double>("ebqe_n");
1544  xt::pyarray<double>& q_r = args.array<double>("q_r");
1545  xt::pyarray<double>& q_porosity = args.array<double>("q_porosity");
1546  int offset_u = args.scalar<int>("offset_u");
1547  int stride_u = args.scalar<int>("stride_u");
1548  xt::pyarray<double>& globalResidual = args.array<double>("globalResidual");
1549  int nExteriorElementBoundaries_global = args.scalar<int>("nExteriorElementBoundaries_global");
1550  xt::pyarray<int>& exteriorElementBoundariesArray = args.array<int>("exteriorElementBoundariesArray");
1551  xt::pyarray<int>& elementBoundaryElementsArray = args.array<int>("elementBoundaryElementsArray");
1552  xt::pyarray<int>& elementBoundaryLocalElementBoundariesArray = args.array<int>("elementBoundaryLocalElementBoundariesArray");
1553  xt::pyarray<double>& H_dof = args.array<double>("H_dof");
1554  gf.useExact=useExact;
1555  gf_nodes.useExact=useExact;
1556  for(int eN=0;eN<nElements_global;eN++)
1557  {
1558  double epsHeaviside;
1559  //loop over quadrature points and compute integrands
1560  //declare local storage for element residual and initialize
1561  register double element_phi[nDOF_trial_element];
1562  for (int i=0;i<nDOF_test_element;i++)
1563  {
1564  register int eN_i=eN*nDOF_test_element+i;
1565  element_phi[i] = phi_dof.data()[phi_l2g.data()[eN_i]];
1566  }//i
1567  double element_nodes[nDOF_mesh_trial_element*3];
1568  for (int i=0;i<nDOF_mesh_trial_element;i++)
1569  {
1570  register int eN_i=eN*nDOF_mesh_trial_element+i;
1571  for(int I=0;I<3;I++)
1572  element_nodes[i*3 + I] = mesh_dof.data()[mesh_l2g.data()[eN_i]*3 + I];
1573  }//i
1574  gf.calculate(element_phi, element_nodes, x_ref.data(),false);
1575  gf_nodes.calculate(element_phi, element_nodes, element_nodes,false);
1576  for (int k=0;k<nQuadraturePoints_element;k++)
1577  {
1578  //compute indeces and declare local storage
1579  register int eN_k = eN*nQuadraturePoints_element+k,
1580  eN_k_nSpace = eN_k*nSpace;
1581  //eN_nDOF_trial_element = eN*nDOF_trial_element;
1582  //register double u=0.0,grad_u[nSpace],r=0.0,dr=0.0;
1583  register double jac[nSpace*nSpace],
1584  jacDet,
1585  jacInv[nSpace*nSpace],
1586  //u_grad_trial[nDOF_trial_element*nSpace],
1587  //u_test_dV[nDOF_trial_element],
1588  //u_grad_test_dV[nDOF_test_element*nSpace],
1589  dV,x,y,z,
1590  G[nSpace*nSpace],G_dd_G,tr_G,h_phi;
1591  //
1592  //compute solution and gradients at quadrature points
1593  //
1594  gf.set_quad(k);
1595  ck.calculateMapping_element(eN,
1596  k,
1597  mesh_dof.data(),
1598  mesh_l2g.data(),
1599  mesh_trial_ref.data(),
1600  mesh_grad_trial_ref.data(),
1601  jac,
1602  jacDet,
1603  jacInv,
1604  x,y,z);
1605  ck.calculateH_element(eN,
1606  k,
1607  nodeDiametersArray.data(),
1608  mesh_l2g.data(),
1609  mesh_trial_ref.data(),
1610  h_phi);
1611  //get the physical integration weight
1612  dV = fabs(jacDet)*dV_ref.data()[k];
1613  ck.calculateG(jacInv,G,G_dd_G,tr_G);
1614 
1615  /* double dir[nSpace]; */
1616  /* double norm = 1.0e-8; */
1617  /* for (int I=0;I<nSpace;I++) */
1618  /* norm += q_normal_phi.data()[eN_k_nSpace+I]*q_normal_phi.data()[eN_k_nSpace+I]; */
1619  /* norm = sqrt(norm); */
1620  /* for (int I=0;I<nSpace;I++) */
1621  /* dir[I] = q_normal_phi.data()[eN_k_nSpace+I]/norm; */
1622 
1623  /* ck.calculateGScale(G,dir,h_phi); */
1624  epsHeaviside=epsFactHeaviside*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter.data()[eN]);
1625  q_H.data()[eN_k] = gf.H(epsHeaviside,q_phi.data()[eN_k]);
1626  }//k
1627  // distribute rhs for mass correction
1628  for (int i=0;i<nDOF_trial_element;i++)
1629  {
1630  gf_nodes.set_quad(i);
1631  int eN_i = eN*nDOF_trial_element + i;
1632  int gi = phi_l2g.data()[eN_i];
1633  epsHeaviside = epsFactHeaviside*nodeDiametersArray.data()[mesh_l2g.data()[eN_i]];//cek hack, only works if isoparametric, but we can fix by including interpolation points
1634  H_dof.data() [gi] = gf_nodes.H(epsHeaviside,phi_dof.data()[gi]);
1635  }
1636  }//elements
1637  }
1638 
1639  void FCTStep(arguments_dict& args)
1640  {
1641  int NNZ = args.scalar<int>("NNZ");
1642  int numDOFs = args.scalar<int>("numDOFs");
1643  xt::pyarray<double>& lumped_mass_matrix = args.array<double>("lumped_mass_matrix");
1644  xt::pyarray<double>& solH = args.array<double>("solH");
1645  xt::pyarray<double>& solL = args.array<double>("solL");
1646  xt::pyarray<double>& limited_solution = args.array<double>("limited_solution");
1647  xt::pyarray<int>& csrRowIndeces_DofLoops = args.array<int>("csrRowIndeces_DofLoops");
1648  xt::pyarray<int>& csrColumnOffsets_DofLoops = args.array<int>("csrColumnOffsets_DofLoops");
1649  xt::pyarray<double>& MassMatrix = args.array<double>("matrix");
1650  Rpos.resize(numDOFs,0.0), Rneg.resize(numDOFs,0.0);
1651  FluxCorrectionMatrix.resize(NNZ,0.0);
1653  // LOOP in DOFs //
1655  int ij=0;
1656  for (int i=0; i<numDOFs; i++)
1657  {
1658  //read some vectors
1659  double solHi = solH.data()[i];
1660  double solLi = solL.data()[i];
1661  double mi = lumped_mass_matrix.data()[i];
1662 
1663  double mini=0., maxi=1.0;
1664  double Pposi=0, Pnegi=0;
1665  // LOOP OVER THE SPARSITY PATTERN (j-LOOP)//
1666  for (int offset=csrRowIndeces_DofLoops.data()[i]; offset<csrRowIndeces_DofLoops.data()[i+1]; offset++)
1667  {
1668  int j = csrColumnOffsets_DofLoops.data()[offset];
1669  // i-th row of flux correction matrix
1670  FluxCorrectionMatrix[ij] = ((i==j ? 1. : 0.)*mi - MassMatrix.data()[ij]) * (solH.data()[j]-solHi);
1671 
1673  // COMPUTE P VECTORS //
1675  Pposi += FluxCorrectionMatrix[ij]*((FluxCorrectionMatrix[ij] > 0) ? 1. : 0.);
1676  Pnegi += FluxCorrectionMatrix[ij]*((FluxCorrectionMatrix[ij] < 0) ? 1. : 0.);
1677 
1678  //update ij
1679  ij+=1;
1680  }
1682  // COMPUTE Q VECTORS //
1684  double Qposi = mi*(maxi-solLi);
1685  double Qnegi = mi*(mini-solLi);
1686 
1688  // COMPUTE R VECTORS //
1690  Rpos[i] = ((Pposi==0) ? 1. : std::min(1.0,Qposi/Pposi));
1691  Rneg[i] = ((Pnegi==0) ? 1. : std::min(1.0,Qnegi/Pnegi));
1692  } // i DOFs
1693 
1695  // COMPUTE LIMITERS //
1697  ij=0;
1698  for (int i=0; i<numDOFs; i++)
1699  {
1700  double ith_Limiter_times_FluxCorrectionMatrix = 0.;
1701  double Rposi = Rpos[i], Rnegi = Rneg[i];
1702  // LOOP OVER THE SPARSITY PATTERN (j-LOOP)//
1703  for (int offset=csrRowIndeces_DofLoops.data()[i]; offset<csrRowIndeces_DofLoops.data()[i+1]; offset++)
1704  {
1705  int j = csrColumnOffsets_DofLoops.data()[offset];
1706  ith_Limiter_times_FluxCorrectionMatrix +=
1707  ((FluxCorrectionMatrix[ij]>0) ? std::min(Rposi,Rneg[j]) : std::min(Rnegi,Rpos[j]))
1708  * FluxCorrectionMatrix[ij];
1709  //ith_Limiter_times_FluxCorrectionMatrix += FluxCorrectionMatrix[ij];
1710  //update ij
1711  ij+=1;
1712  }
1713  limited_solution.data()[i] = fmax(0.0,solL.data()[i] + 1./lumped_mass_matrix.data()[i]*ith_Limiter_times_FluxCorrectionMatrix);
1714  }
1715  }
1716 
1717  // mql. copied from calculateElementJacobian. NOTE: there are some not necessary computations!!!
1718  inline void calculateElementMassMatrix(//element
1719  double* mesh_trial_ref,
1720  double* mesh_grad_trial_ref,
1721  double* mesh_dof,
1722  int* mesh_l2g,
1723  double* dV_ref,
1724  double* u_trial_ref,
1725  double* u_grad_trial_ref,
1726  double* u_test_ref,
1727  double* u_grad_test_ref,
1728  //element boundary
1729  double* mesh_trial_trace_ref,
1730  double* mesh_grad_trial_trace_ref,
1731  double* dS_ref,
1732  double* u_trial_trace_ref,
1733  double* u_grad_trial_trace_ref,
1734  double* u_test_trace_ref,
1735  double* u_grad_test_trace_ref,
1736  double* normal_ref,
1737  double* boundaryJac_ref,
1738  //physics
1739  int nElements_global,
1740  double useMetrics,
1741  double epsFactHeaviside,
1742  double epsFactDirac,
1743  double epsFactDiffusion,
1744  int* u_l2g,
1745  double* elementDiameter,
1746  double* nodeDiametersArray,
1747  double* u_dof,
1748  // double* u_trial,
1749  // double* u_grad_trial,
1750  // double* u_test_dV,
1751  // double* u_grad_test_dV,
1752  double* q_phi,
1753  double* q_normal_phi,
1754  double* q_H,
1755  double* q_porosity,
1756  double* elementMassMatrix,
1757  double* elementLumpedMassMatrix,
1758  double* element_u,
1759  int eN)
1760  {
1761  for (int i=0;i<nDOF_test_element;i++)
1762  {
1763  elementLumpedMassMatrix[i] = 0.0;
1764  for (int j=0;j<nDOF_trial_element;j++)
1765  {
1766  elementMassMatrix[i*nDOF_trial_element+j]=0.0;
1767  }
1768  }
1769  double epsHeaviside,epsDirac,epsDiffusion;
1770  for (int k=0;k<nQuadraturePoints_element;k++)
1771  {
1772  int eN_k = eN*nQuadraturePoints_element+k, //index to a scalar at a quadrature point
1773  eN_k_nSpace = eN_k*nSpace;
1774  //eN_nDOF_trial_element = eN*nDOF_trial_element; //index to a vector at a quadrature point
1775 
1776  //declare local storage
1777  register double u=0.0,
1778  grad_u[nSpace],
1779  r=0.0,dr=0.0,
1780  jac[nSpace*nSpace],
1781  jacDet,
1782  jacInv[nSpace*nSpace],
1783  u_grad_trial[nDOF_trial_element*nSpace],
1784  dV,
1785  u_test_dV[nDOF_test_element],
1786  u_grad_test_dV[nDOF_test_element*nSpace],
1787  x,y,z,
1788  G[nSpace*nSpace],G_dd_G,tr_G,h_phi;
1789  //
1790  //calculate solution and gradients at quadrature points
1791  //
1792  ck.calculateMapping_element(eN,
1793  k,
1794  mesh_dof,
1795  mesh_l2g,
1796  mesh_trial_ref,
1797  mesh_grad_trial_ref,
1798  jac,
1799  jacDet,
1800  jacInv,
1801  x,y,z);
1802  ck.calculateH_element(eN,
1803  k,
1804  nodeDiametersArray,
1805  mesh_l2g,
1806  mesh_trial_ref,
1807  h_phi);
1808  //get the physical integration weight
1809  dV = fabs(jacDet)*dV_ref[k];
1810  ck.calculateG(jacInv,G,G_dd_G,tr_G);
1811 
1812  /* double dir[nSpace]; */
1813  /* double norm = 1.0e-8; */
1814  /* for (int I=0;I<nSpace;I++) */
1815  /* norm += q_normal_phi[eN_k_nSpace+I]*q_normal_phi[eN_k_nSpace+I]; */
1816  /* norm = sqrt(norm); */
1817  /* for (int I=0;I<nSpace;I++) */
1818  /* dir[I] = q_normal_phi[eN_k_nSpace+I]/norm; */
1819  /* ck.calculateGScale(G,dir,h_phi); */
1820 
1821 
1822  //get the trial function gradients
1823  ck.gradTrialFromRef(&u_grad_trial_ref[k*nDOF_trial_element*nSpace],jacInv,u_grad_trial);
1824  //get the solution
1825  ck.valFromElementDOF(element_u,&u_trial_ref[k*nDOF_trial_element],u);
1826  //get the solution gradients
1827  ck.gradFromElementDOF(element_u,u_grad_trial,grad_u);
1828  //precalculate test function products with integration weights
1829  for (int j=0;j<nDOF_trial_element;j++)
1830  {
1831  u_test_dV[j] = u_test_ref[k*nDOF_trial_element+j]*dV;
1832  for (int I=0;I<nSpace;I++)
1833  {
1834  u_grad_test_dV[j*nSpace+I] = u_grad_trial[j*nSpace+I]*dV;//cek warning won't work for Petrov-Galerkin
1835  }
1836  }
1837  //
1838  //calculate pde coefficients and derivatives at quadrature points
1839  //
1840  epsHeaviside=epsFactHeaviside*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter[eN]);
1841  epsDirac =epsFactDirac* (useMetrics*h_phi+(1.0-useMetrics)*elementDiameter[eN]);
1842  epsDiffusion=epsFactDiffusion*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter[eN]);
1843  // *(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter[eN]);
1844  evaluateCoefficients(epsHeaviside,
1845  epsDirac,
1846  q_phi[eN_k],
1847  q_H[eN_k],
1848  u,
1849  q_porosity[eN_k],
1850  r,
1851  dr);
1852  for(int i=0;i<nDOF_test_element;i++)
1853  {
1854  //int eN_k_i=eN_k*nDOF_test_element+i;
1855  //int eN_k_i_nSpace=eN_k_i*nSpace;
1856  elementLumpedMassMatrix[i] += u_test_dV[i];
1857  for(int j=0;j<nDOF_trial_element;j++)
1858  {
1859  elementMassMatrix[i*nDOF_trial_element+j] += u_trial_ref[k*nDOF_trial_element+j]*u_test_dV[i];
1860  }//j
1861  }//i
1862  }//k
1863  }
1864 
1865  void calculateMassMatrix(arguments_dict& args)
1866  {
1867  xt::pyarray<double>& mesh_trial_ref = args.array<double>("mesh_trial_ref");
1868  xt::pyarray<double>& mesh_grad_trial_ref = args.array<double>("mesh_grad_trial_ref");
1869  xt::pyarray<double>& mesh_dof = args.array<double>("mesh_dof");
1870  xt::pyarray<int>& mesh_l2g = args.array<int>("mesh_l2g");
1871  xt::pyarray<double>& dV_ref = args.array<double>("dV_ref");
1872  xt::pyarray<double>& u_trial_ref = args.array<double>("u_trial_ref");
1873  xt::pyarray<double>& u_grad_trial_ref = args.array<double>("u_grad_trial_ref");
1874  xt::pyarray<double>& u_test_ref = args.array<double>("u_test_ref");
1875  xt::pyarray<double>& u_grad_test_ref = args.array<double>("u_grad_test_ref");
1876  xt::pyarray<double>& mesh_trial_trace_ref = args.array<double>("mesh_trial_trace_ref");
1877  xt::pyarray<double>& mesh_grad_trial_trace_ref = args.array<double>("mesh_grad_trial_trace_ref");
1878  xt::pyarray<double>& dS_ref = args.array<double>("dS_ref");
1879  xt::pyarray<double>& u_trial_trace_ref = args.array<double>("u_trial_trace_ref");
1880  xt::pyarray<double>& u_grad_trial_trace_ref = args.array<double>("u_grad_trial_trace_ref");
1881  xt::pyarray<double>& u_test_trace_ref = args.array<double>("u_test_trace_ref");
1882  xt::pyarray<double>& u_grad_test_trace_ref = args.array<double>("u_grad_test_trace_ref");
1883  xt::pyarray<double>& normal_ref = args.array<double>("normal_ref");
1884  xt::pyarray<double>& boundaryJac_ref = args.array<double>("boundaryJac_ref");
1885  int nElements_global = args.scalar<int>("nElements_global");
1886  double useMetrics = args.scalar<double>("useMetrics");
1887  double epsFactHeaviside = args.scalar<double>("epsFactHeaviside");
1888  double epsFactDirac = args.scalar<double>("epsFactDirac");
1889  double epsFactDiffusion = args.scalar<double>("epsFactDiffusion");
1890  xt::pyarray<int>& u_l2g = args.array<int>("u_l2g");
1891  xt::pyarray<double>& elementDiameter = args.array<double>("elementDiameter");
1892  xt::pyarray<double>& nodeDiametersArray = args.array<double>("nodeDiametersArray");
1893  xt::pyarray<double>& u_dof = args.array<double>("u_dof");
1894  xt::pyarray<double>& q_phi = args.array<double>("q_phi");
1895  xt::pyarray<double>& q_normal_phi = args.array<double>("q_normal_phi");
1896  xt::pyarray<double>& q_H = args.array<double>("q_H");
1897  xt::pyarray<double>& q_porosity = args.array<double>("q_porosity");
1898  xt::pyarray<int>& csrRowIndeces_u_u = args.array<int>("csrRowIndeces_u_u");
1899  xt::pyarray<int>& csrColumnOffsets_u_u = args.array<int>("csrColumnOffsets_u_u");
1900  xt::pyarray<double>& globalMassMatrix = args.array<double>("globalMassMatrix");
1901  xt::pyarray<double>& globalLumpedMassMatrix = args.array<double>("globalLumpedMassMatrix");
1902  //
1903  //loop over elements to compute volume integrals and load them into the element Jacobians and global Jacobian
1904  //
1905  for(int eN=0;eN<nElements_global;eN++)
1906  {
1907  register double elementMassMatrix[nDOF_test_element*nDOF_trial_element],element_u[nDOF_trial_element], elementLumpedMassMatrix[nDOF_trial_element];
1908  for (int j=0;j<nDOF_trial_element;j++)
1909  {
1910  register int eN_j = eN*nDOF_trial_element+j;
1911  element_u[j] = u_dof.data()[u_l2g.data()[eN_j]];
1912  }
1913  calculateElementMassMatrix(mesh_trial_ref.data(),
1914  mesh_grad_trial_ref.data(),
1915  mesh_dof.data(),
1916  mesh_l2g.data(),
1917  dV_ref.data(),
1918  u_trial_ref.data(),
1919  u_grad_trial_ref.data(),
1920  u_test_ref.data(),
1921  u_grad_test_ref.data(),
1922  mesh_trial_trace_ref.data(),
1923  mesh_grad_trial_trace_ref.data(),
1924  dS_ref.data(),
1925  u_trial_trace_ref.data(),
1926  u_grad_trial_trace_ref.data(),
1927  u_test_trace_ref.data(),
1928  u_grad_test_trace_ref.data(),
1929  normal_ref.data(),
1930  boundaryJac_ref.data(),
1931  nElements_global,
1932  useMetrics,
1933  epsFactHeaviside,
1934  epsFactDirac,
1935  epsFactDiffusion,
1936  u_l2g.data(),
1937  elementDiameter.data(),
1938  nodeDiametersArray.data(),
1939  u_dof.data(),
1940  q_phi.data(),
1941  q_normal_phi.data(),
1942  q_H.data(),
1943  q_porosity.data(),
1944  elementMassMatrix,
1945  elementLumpedMassMatrix,
1946  element_u,
1947  eN);
1948  //
1949  //load into element Jacobian into global Jacobian
1950  //
1951  for (int i=0;i<nDOF_test_element;i++)
1952  {
1953  int eN_i = eN*nDOF_test_element+i;
1954  int gi = u_l2g.data()[eN_i];
1955  globalLumpedMassMatrix.data()[gi] += elementLumpedMassMatrix[i];
1956  for (int j=0;j<nDOF_trial_element;j++)
1957  {
1958  int eN_i_j = eN_i*nDOF_trial_element+j;
1959  globalMassMatrix.data()[csrRowIndeces_u_u.data()[eN_i] + csrColumnOffsets_u_u.data()[eN_i_j]] +=
1960  elementMassMatrix[i*nDOF_trial_element+j];
1961  }//j
1962  }//i
1963  }//elements
1964  }//calculate mass matrix
1965 
1966  void setMassQuadratureEdgeBasedStabilizationMethods(arguments_dict& args,
1967  bool useExact)
1968  {
1969  xt::pyarray<double>& mesh_trial_ref = args.array<double>("mesh_trial_ref");
1970  xt::pyarray<double>& mesh_grad_trial_ref = args.array<double>("mesh_grad_trial_ref");
1971  xt::pyarray<double>& mesh_dof = args.array<double>("mesh_dof");
1972  xt::pyarray<int>& mesh_l2g = args.array<int>("mesh_l2g");
1973  xt::pyarray<double>& x_ref = args.array<double>("x_ref");
1974  xt::pyarray<double>& dV_ref = args.array<double>("dV_ref");
1975  xt::pyarray<double>& u_trial_ref = args.array<double>("u_trial_ref");
1976  xt::pyarray<double>& u_grad_trial_ref = args.array<double>("u_grad_trial_ref");
1977  xt::pyarray<double>& u_test_ref = args.array<double>("u_test_ref");
1978  xt::pyarray<double>& u_grad_test_ref = args.array<double>("u_grad_test_ref");
1979  xt::pyarray<double>& mesh_trial_trace_ref = args.array<double>("mesh_trial_trace_ref");
1980  xt::pyarray<double>& mesh_grad_trial_trace_ref = args.array<double>("mesh_grad_trial_trace_ref");
1981  xt::pyarray<double>& dS_ref = args.array<double>("dS_ref");
1982  xt::pyarray<double>& u_trial_trace_ref = args.array<double>("u_trial_trace_ref");
1983  xt::pyarray<double>& u_grad_trial_trace_ref = args.array<double>("u_grad_trial_trace_ref");
1984  xt::pyarray<double>& u_test_trace_ref = args.array<double>("u_test_trace_ref");
1985  xt::pyarray<double>& u_grad_test_trace_ref = args.array<double>("u_grad_test_trace_ref");
1986  xt::pyarray<double>& normal_ref = args.array<double>("normal_ref");
1987  xt::pyarray<double>& boundaryJac_ref = args.array<double>("boundaryJac_ref");
1988  int nElements_global = args.scalar<int>("nElements_global");
1989  double useMetrics = args.scalar<double>("useMetrics");
1990  double epsFactHeaviside = args.scalar<double>("epsFactHeaviside");
1991  double epsFactDirac = args.scalar<double>("epsFactDirac");
1992  double epsFactDiffusion = args.scalar<double>("epsFactDiffusion");
1993  xt::pyarray<int>& phi_l2g = args.array<int>("phi_l2g");
1994  xt::pyarray<double>& elementDiameter = args.array<double>("elementDiameter");
1995  xt::pyarray<double>& nodeDiametersArray = args.array<double>("nodeDiametersArray");
1996  xt::pyarray<double>& phi_dof = args.array<double>("phi_dof");
1997  xt::pyarray<double>& q_phi = args.array<double>("q_phi");
1998  xt::pyarray<double>& q_normal_phi = args.array<double>("q_normal_phi");
1999  xt::pyarray<double>& ebqe_phi = args.array<double>("ebqe_phi");
2000  xt::pyarray<double>& ebqe_normal_phi = args.array<double>("ebqe_normal_phi");
2001  xt::pyarray<double>& q_H = args.array<double>("q_H");
2002  xt::pyarray<double>& q_u = args.array<double>("q_u");
2003  xt::pyarray<double>& q_n = args.array<double>("q_n");
2004  xt::pyarray<double>& ebqe_u = args.array<double>("ebqe_u");
2005  xt::pyarray<double>& ebqe_n = args.array<double>("ebqe_n");
2006  xt::pyarray<double>& q_r = args.array<double>("q_r");
2007  xt::pyarray<double>& q_porosity = args.array<double>("q_porosity");
2008  int offset_u = args.scalar<int>("offset_u");
2009  int stride_u = args.scalar<int>("stride_u");
2010  xt::pyarray<double>& globalResidual = args.array<double>("globalResidual");
2011  int nExteriorElementBoundaries_global = args.scalar<int>("nExteriorElementBoundaries_global");
2012  xt::pyarray<int>& exteriorElementBoundariesArray = args.array<int>("exteriorElementBoundariesArray");
2013  xt::pyarray<int>& elementBoundaryElementsArray = args.array<int>("elementBoundaryElementsArray");
2014  xt::pyarray<int>& elementBoundaryLocalElementBoundariesArray = args.array<int>("elementBoundaryLocalElementBoundariesArray");
2015  xt::pyarray<double>& rhs_mass_correction = args.array<double>("rhs_mass_correction");
2016  xt::pyarray<double>& lumped_L2p_vof_mass_correction = args.array<double>("lumped_L2p_vof_mass_correction");
2017  xt::pyarray<double>& lumped_mass_matrix = args.array<double>("lumped_mass_matrix");
2018  int numDOFs = args.scalar<int>("numDOFs");
2019  gf.useExact=useExact;
2020  for(int eN=0;eN<nElements_global;eN++)
2021  {
2022  register double element_rhs_mass_correction[nDOF_test_element];
2023  for (int i=0;i<nDOF_test_element;i++)
2024  element_rhs_mass_correction[i] = 0.;
2025  double epsHeaviside;
2026  //loop over quadrature points and compute integrands
2027  double element_phi[nDOF_trial_element];
2028  for (int i=0;i<nDOF_test_element;i++)
2029  {
2030  register int eN_i=eN*nDOF_test_element+i;
2031  element_phi[i] = phi_dof.data()[phi_l2g.data()[eN_i]];
2032  }//i
2033  double element_nodes[nDOF_mesh_trial_element*3];
2034  for (int i=0;i<nDOF_mesh_trial_element;i++)
2035  {
2036  register int eN_i=eN*nDOF_mesh_trial_element+i;
2037  for(int I=0;I<3;I++)
2038  element_nodes[i*3 + I] = mesh_dof.data()[mesh_l2g.data()[eN_i]*3 + I];
2039  }//i
2040  gf.calculate(element_phi, element_nodes, x_ref.data(),false);
2041  for (int k=0;k<nQuadraturePoints_element;k++)
2042  {
2043  //compute indeces and declare local storage
2044  register int eN_k = eN*nQuadraturePoints_element+k,
2045  eN_k_nSpace = eN_k*nSpace;
2046  //eN_nDOF_trial_element = eN*nDOF_trial_element;
2047  //register double u=0.0,grad_u[nSpace],r=0.0,dr=0.0;
2048  register double jac[nSpace*nSpace],
2049  jacDet,
2050  jacInv[nSpace*nSpace],
2051  //u_grad_trial[nDOF_trial_element*nSpace],
2052  //u_test_dV[nDOF_trial_element],
2053  //u_grad_test_dV[nDOF_test_element*nSpace],
2054  dV,x,y,z,
2055  u_test_dV[nDOF_test_element],
2056  G[nSpace*nSpace],G_dd_G,tr_G,h_phi;
2057  gf.set_quad(k);
2058  //
2059  ck.calculateMapping_element(eN,
2060  k,
2061  mesh_dof.data(),
2062  mesh_l2g.data(),
2063  mesh_trial_ref.data(),
2064  mesh_grad_trial_ref.data(),
2065  jac,
2066  jacDet,
2067  jacInv,
2068  x,y,z);
2069  ck.calculateH_element(eN,
2070  k,
2071  nodeDiametersArray.data(),
2072  mesh_l2g.data(),
2073  mesh_trial_ref.data(),
2074  h_phi);
2075  //get the physical integration weight
2076  dV = fabs(jacDet)*dV_ref.data()[k];
2077  ck.calculateG(jacInv,G,G_dd_G,tr_G);
2078 
2079  // precalculate test function times integration weight
2080  for (int j=0;j<nDOF_trial_element;j++)
2081  u_test_dV[j] = u_test_ref.data()[k*nDOF_trial_element+j]*dV;
2082 
2083  /* double dir[nSpace]; */
2084  /* double norm = 1.0e-8; */
2085  /* for (int I=0;I<nSpace;I++) */
2086  /* norm += q_normal_phi.data()[eN_k_nSpace+I]*q_normal_phi.data()[eN_k_nSpace+I]; */
2087  /* norm = sqrt(norm); */
2088  /* for (int I=0;I<nSpace;I++) */
2089  /* dir[I] = q_normal_phi.data()[eN_k_nSpace+I]/norm; */
2090 
2091  /* ck.calculateGScale(G,dir,h_phi); */
2092  epsHeaviside=epsFactHeaviside*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter.data()[eN]);
2093  q_H.data()[eN_k] = gf.H(epsHeaviside,q_phi.data()[eN_k]);
2094 
2095  for (int i=0;i<nDOF_trial_element;i++)
2096  element_rhs_mass_correction [i] += q_porosity.data()[eN_k]*q_H.data()[eN_k]*u_test_dV[i];
2097  }//k
2098  // distribute rhs for mass correction
2099  for (int i=0;i<nDOF_trial_element;i++)
2100  {
2101  int eN_i = eN*nDOF_trial_element + i;
2102  int gi = phi_l2g.data()[eN_i];
2103  rhs_mass_correction.data()[gi] += element_rhs_mass_correction[i];
2104  }
2105  }//elements
2106  // COMPUTE LUMPED L2 PROYJECTION
2107  for (int i=0; i<numDOFs; i++)
2108  {
2109  double mi = lumped_mass_matrix.data()[i];
2110  lumped_L2p_vof_mass_correction.data()[i] = 1./mi*rhs_mass_correction.data()[i];
2111  }
2112  }
2113  };//MCorr
2114 
2115  inline MCorr_base* newMCorr(int nSpaceIn,
2116  int nQuadraturePoints_elementIn,
2117  int nDOF_mesh_trial_elementIn,
2118  int nDOF_trial_elementIn,
2119  int nDOF_test_elementIn,
2120  int nQuadraturePoints_elementBoundaryIn,
2121  int CompKernelFlag)
2122  {
2123  if (nSpaceIn == 2)
2124  return proteus::chooseAndAllocateDiscretization2D<MCorr_base,MCorr,CompKernel>(nSpaceIn,
2125  nQuadraturePoints_elementIn,
2126  nDOF_mesh_trial_elementIn,
2127  nDOF_trial_elementIn,
2128  nDOF_test_elementIn,
2129  nQuadraturePoints_elementBoundaryIn,
2130  CompKernelFlag);
2131  else
2132  return proteus::chooseAndAllocateDiscretization<MCorr_base,MCorr,CompKernel>(nSpaceIn,
2133  nQuadraturePoints_elementIn,
2134  nDOF_mesh_trial_elementIn,
2135  nDOF_trial_elementIn,
2136  nDOF_test_elementIn,
2137  nQuadraturePoints_elementBoundaryIn,
2138  CompKernelFlag);
2139  }
2140 }//proteus
2141 #endif
proteus::MCorr::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_porosity, int offset_u, int stride_u, double *elementResidual_u, int nExteriorElementBoundaries_global, int *exteriorElementBoundariesArray, int *elementBoundaryElementsArray, int *elementBoundaryLocalElementBoundariesArray, double *element_u, int eN)
Definition: MCorr.h:72
proteus::MCorr::elementConstantSolve
void elementConstantSolve(arguments_dict &args)
Definition: MCorr.h:994
proteus::MCorr::calculateResidual
void calculateResidual(arguments_dict &args, bool useExact)
Definition: MCorr.h:239
proteus::MCorr::globalConstantRJ
std::tuple< double, double > globalConstantRJ(arguments_dict &args)
Definition: MCorr.h:1222
proteus::MCorr_base::globalConstantRJ
virtual std::tuple< double, double > globalConstantRJ(arguments_dict &args)=0
proteus::MCorr_base::elementSolve
virtual void elementSolve(arguments_dict &args)=0
proteus::MCorr::setMassQuadratureEdgeBasedStabilizationMethods
void setMassQuadratureEdgeBasedStabilizationMethods(arguments_dict &args, bool useExact)
Definition: MCorr.h:1966
proteus::MCorr::setMassQuadrature
void setMassQuadrature(arguments_dict &args, bool useExact)
Definition: MCorr.h:1504
proteus::MCorr_base::calculateResidual
virtual void calculateResidual(arguments_dict &args, bool useExact)=0
proteus::MCorr_base::calculateMassMatrix
virtual void calculateMassMatrix(arguments_dict &args)=0
proteus::MCorr_base::setMassQuadrature
virtual void setMassQuadrature(arguments_dict &args, bool useExact)=0
dgesc2_
int dgesc2_(int *n, double *a, int *lda, double *rhs, int *ipiv, int *jpiv, double *scale)
proteus::MCorr_base
Definition: MCorr.h:24
proteus::MCorr_base::FluxCorrectionMatrix
std::valarray< double > FluxCorrectionMatrix
Definition: MCorr.h:27
proteus::MCorr::calculateMassMatrix
void calculateMassMatrix(arguments_dict &args)
Definition: MCorr.h:1865
proteus::MCorr::FCTStep
void FCTStep(arguments_dict &args)
Definition: MCorr.h:1639
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::MCorr::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_porosity, double *elementJacobian_u_u, double *element_u, int eN)
Definition: MCorr.h:467
proteus::MCorr::gf
GeneralizedFunctions< nSpace, 2, nQuadraturePoints_element, nQuadraturePoints_elementBoundary > gf
Definition: MCorr.h:53
H
Double H
Definition: Headers.h:65
min
#define min(a, b)
Definition: jf.h:71
proteus::MCorr
Definition: MCorr.h:50
proteus::MCorr_base::elementConstantSolve
virtual void elementConstantSolve(arguments_dict &args)=0
proteus::MCorr::ck
CompKernelType ck
Definition: MCorr.h:52
proteus::MCorr::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: MCorr.h:59
equivalent_polynomials.h
z
Double * z
Definition: Headers.h:49
u
Double u
Definition: Headers.h:89
proteus::MCorr_base::setMassQuadratureEdgeBasedStabilizationMethods
virtual void setMassQuadratureEdgeBasedStabilizationMethods(arguments_dict &args, bool useExact)=0
proteus::MCorr_base::FCTStep
virtual void FCTStep(arguments_dict &args)=0
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::MCorr_base::~MCorr_base
virtual ~MCorr_base()
Definition: MCorr.h:28
proteus::MCorr_base::calculateJacobian
virtual void calculateJacobian(arguments_dict &args, bool useExact)=0
dgetc2_
int dgetc2_(int *n, double *a, int *lda, int *ipiv, int *jpiv, int *info)
equivalent_polynomials::GeneralizedFunctions_mix
Definition: equivalent_polynomials.h:767
proteus::MCorr::elementSolve
void elementSolve(arguments_dict &args)
Definition: MCorr.h:724
proteus::MCorr_base::Rneg
std::valarray< double > Rneg
Definition: MCorr.h:26
proteus
Definition: ADR.h:17
proteus::GeneralizedFunctions
equivalent_polynomials::GeneralizedFunctions_mix< nSpace, nP, nQ, nEBQ > GeneralizedFunctions
Definition: ADR.h:19
proteus::MCorr::calculateElementMassMatrix
void calculateElementMassMatrix(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_porosity, double *elementMassMatrix, double *elementLumpedMassMatrix, double *element_u, int eN)
Definition: MCorr.h:1718
proteus::MCorr_base::Rpos
std::valarray< double > Rpos
Definition: MCorr.h:26
r
Double r
Definition: Headers.h:83
proteus::MCorr::gf_nodes
GeneralizedFunctions< nSpace, 2, nDOF_trial_element, nQuadraturePoints_elementBoundary > gf_nodes
Definition: MCorr.h:54
proteus::arguments_dict
Definition: ArgumentsDict.h:70
proteus::MCorr_base::calculateMass
virtual double calculateMass(arguments_dict &args, bool useExact)=0
proteus::MCorr::MCorr
MCorr()
Definition: MCorr.h:55
proteus::MCorr::calculateMass
double calculateMass(arguments_dict &args, bool useExact)
Definition: MCorr.h:1380
ArgumentsDict.h
proteus::newMCorr
MCorr_base * newMCorr(int nSpaceIn, int nQuadraturePoints_elementIn, int nDOF_mesh_trial_elementIn, int nDOF_trial_elementIn, int nDOF_test_elementIn, int nQuadraturePoints_elementBoundaryIn, int CompKernelFlag)
Definition: MCorr.h:2115
proteus::MCorr::calculateJacobian
void calculateJacobian(arguments_dict &args, bool useExact)
Definition: MCorr.h:615
Generated on Fri Jul 1 2022 11:17:34 for proteus by doxygen 1.8.20