proteus 1.9.0
C/C++/Fortran libraries
Loading...
Searching...
No Matches
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
14namespace 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 int eN_k = eN*nQuadraturePoints_element+k,
162 eN_k_nSpace = eN_k*nSpace;
163 //eN_nDOF_trial_element = eN*nDOF_trial_element;
164 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 //int eN_k_i=eN_k*nDOF_test_element+i;
241 //int eN_k_i_nSpace = eN_k_i*nSpace;
242 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 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 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 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 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 //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 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 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 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 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 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 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 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 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 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 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 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 int eN_k = eN*nQuadraturePoints_element+k,
1584 eN_k_nSpace = eN_k*nSpace;
1585 //eN_nDOF_trial_element = eN*nDOF_trial_element;
1586 //double u=0.0,grad_u[nSpace],r=0.0,dr=0.0;
1587 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 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 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 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
HI
#define HI
Definition Headers.h:4
r
Double r
Definition Headers.h:83
H
Double H
Definition Headers.h:65
u
Double u
Definition Headers.h:89
z
Double * z
Definition Headers.h:49
FAST_ASSEMBLY
#define FAST_ASSEMBLY
Definition MCorr3P.h:12
proteus::cppMCorr3P_base
Definition MCorr3P.h:17
proteus::cppMCorr3P_base::~cppMCorr3P_base
virtual ~cppMCorr3P_base()
Definition MCorr3P.h:19
proteus::cppMCorr3P_base::elementSolve
virtual void elementSolve(arguments_dict &args)=0
proteus::cppMCorr3P_base::setMassQuadrature
virtual void setMassQuadrature(arguments_dict &args)=0
proteus::cppMCorr3P_base::elementConstantSolve
virtual void elementConstantSolve(arguments_dict &args)=0
proteus::cppMCorr3P_base::calculateStiffnessMatrix
virtual void calculateStiffnessMatrix(arguments_dict &args)=0
proteus::cppMCorr3P_base::globalConstantRJ
virtual std::pair< double, double > globalConstantRJ(arguments_dict &args)=0
proteus::cppMCorr3P_base::calculateResidual
virtual void calculateResidual(arguments_dict &args)=0
proteus::cppMCorr3P_base::calculateMass
virtual double calculateMass(arguments_dict &args)=0
proteus::cppMCorr3P_base::calculateJacobian
virtual void calculateJacobian(arguments_dict &args)=0
proteus::cppMCorr3P::elementConstantSolve
void elementConstantSolve(arguments_dict &args)
Definition MCorr3P.h:1032
proteus::cppMCorr3P::smoothedHeaviside_integral
double smoothedHeaviside_integral(double eps, double phi)
Definition MCorr3P.h:57
proteus::cppMCorr3P::calculateStiffnessMatrix
void calculateStiffnessMatrix(arguments_dict &args)
Definition MCorr3P.h:1645
proteus::cppMCorr3P::calculateResidual
void calculateResidual(arguments_dict &args)
Definition MCorr3P.h:267
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
proteus::cppMCorr3P::smoothedHeaviside
double smoothedHeaviside(double eps, double phi)
Definition MCorr3P.h:43
proteus::cppMCorr3P::smoothedDirac
double smoothedDirac(double eps, double phi)
Definition MCorr3P.h:75
proteus::cppMCorr3P::setMassQuadrature
void setMassQuadrature(arguments_dict &args)
Definition MCorr3P.h:1527
proteus::cppMCorr3P::ck
CompKernelType ck
Definition MCorr3P.h:40
proteus::cppMCorr3P::calculateMass
double calculateMass(arguments_dict &args)
Definition MCorr3P.h:1424
proteus::cppMCorr3P::globalConstantRJ
std::pair< double, double > globalConstantRJ(arguments_dict &args)
Definition MCorr3P.h:1264
proteus::cppMCorr3P::elementSolve
void elementSolve(arguments_dict &args)
Definition MCorr3P.h:759
proteus::cppMCorr3P::calculateJacobian
void calculateJacobian(arguments_dict &args)
Definition MCorr3P.h:641
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
Definition ADR.h:19
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::chooseAndAllocateDiscretization
Model_Base * chooseAndAllocateDiscretization(int nSpaceIn, int nQuadraturePoints_elementIn, int nDOF_mesh_trial_elementIn, int nDOF_trial_elementIn, int nDOF_test_elementIn, int nDOF_v_trial_elementIn, int nDOF_v_test_elementIn, int nQuadraturePoints_elementBoundaryIn, int CompKernelFlag)
Definition MixedModelFactory.h:48
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::chooseAndAllocateDiscretization2D
Model_Base * chooseAndAllocateDiscretization2D(int nSpaceIn, int nQuadraturePoints_elementIn, int nDOF_mesh_trial_elementIn, int nDOF_trial_elementIn, int nDOF_test_elementIn, int nDOF_v_trial_elementIn, int nDOF_v_test_elementIn, int nQuadraturePoints_elementBoundaryIn, int CompKernelFlag)
Definition MixedModelFactory.h:309
dgesc2_
int dgesc2_(int *n, double *a, int *lda, double *rhs, int *ipiv, int *jpiv, double *scale)
dgetc2_
int dgetc2_(int *n, double *a, int *lda, int *ipiv, int *jpiv, int *info)
proteus::arguments_dict
Definition ArgumentsDict.h:70
proteus::arguments_dict::scalar
T & scalar(const std::string &key)
proteus::arguments_dict::array
xt::pyarray< T > & array(const std::string &key)
Generated on for proteus by doxygen 1.18.0