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