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