proteus 1.9.0
C/C++/Fortran libraries
Loading...
Searching...
No Matches
PresInc.h
Go to the documentation of this file.
1#ifndef PRESINC_H
2#define PRESINC_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 cppPresInc_base
15 {
16 public:
17 virtual ~cppPresInc_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 cppPresInc : public cppPresInc_base
30 {
31 public:
32 const int nDOF_test_X_trial_element;
33 CompKernelType ck;
34 cppPresInc():
35 nDOF_test_X_trial_element(nDOF_test_element*nDOF_trial_element),
36 ck()
37 {}
38 inline
39 void evaluateCoefficients(const double& alphaBDF,
40 const double vf[nSpace],
41 const double vs[nSpace],
42 const double& vos,
43 const double& rhos_min,
44 const double& rhof_min,
45 double f[nSpace],
46 double& a)
47 {
48 for (int I=0;I<nSpace;I++)
49 f[I] = (1.0-vos)*vf[I] + vos*vs[I];
50 //a = 1.0/( vos*rhos_min*alphaBDF) + 1.0/( (1.0-vos)*rhof_min*alphaBDF );
51 a = 1.0/( vos*rhos_min*alphaBDF + (1.0-vos)*rhof_min*alphaBDF );
52 //a = (1.0-vos)/(rhof_min*alphaBDF) + vos*(rhos_min/alphaBDF);
53 //a = (1.0)/(rhof_min*alphaBDF) ;
54 //a = (1.0-vos)/(rhof_min*alphaBDF) + (vos)/(rhos_min*alphaBDF);
55 }
56
57 inline
58 void exteriorNumericalAdvectiveFlux(const int& isFluxBoundary,
59 const double& bc_flux,
60 const double n[nSpace],
61 const double f[nSpace],
62 double& flux)
63 {
64 if (isFluxBoundary == 1)
65 flux = bc_flux;
66 else
67 {
68 flux = 0.0;
69 for (int I=0; I < nSpace; I++)
70 flux += n[I]*f[I];
71 }
72 }
73
74 inline
75 void exteriorNumericalDiffusiveFlux(const int& isDOFBoundary,
76 const int& isFluxBoundary,
77 const double n[nSpace],
78 const double& a,
79 const double grad_potential[nSpace],
80 const double& u,
81 const double& bc_u,
82 const double& bc_flux,
83 const double& penalty,
84 double& flux)
85 {
86 if(isFluxBoundary == 1)
87 {
88 flux = bc_flux;
89 }
90 else if(isDOFBoundary == 1)
91 {
92 flux = 0.0;
93 for(int I=0;I<nSpace;I++)
94 flux-= a*grad_potential[I]*n[I];
95 flux += a*penalty*(u-bc_u);
96 }
97 else
98 {
99 std::cerr<<"PresInc.h: warning, diffusion term with no boundary condition set, setting diffusive flux to 0.0"<<std::endl;
100 flux = 0.0;
101 }
102 }
103
104 inline
105 double ExteriorNumericalDiffusiveFluxJacobian(const int& isDOFBoundary,
106 const int& isFluxBoundary,
107 const double n[nSpace],
108 const double& a,
109 const double& v,
110 const double grad_v[nSpace],
111 const double& penalty)
112 {
113 double tmp=0.0;
114 if(isFluxBoundary==0 && isDOFBoundary==1)
115 {
116 for(int I=0;I<nSpace;I++)
117 tmp -= a*grad_v[I]*n[I];
118 tmp +=a*penalty*v;
119 }
120 return tmp;
121 }
122
123 inline void calculateElementResidual(//element
124 double* mesh_trial_ref,
125 double* mesh_grad_trial_ref,
126 double* mesh_dof,
127 int* mesh_l2g,
128 double* dV_ref,
129 double* u_trial_ref,
130 double* u_grad_trial_ref,
131 double* u_test_ref,
132 double* u_grad_test_ref,
133 //element boundary
134 double* mesh_trial_trace_ref,
135 double* mesh_grad_trial_trace_ref,
136 double* dS_ref,
137 double* u_trial_trace_ref,
138 double* u_grad_trial_trace_ref,
139 double* u_test_trace_ref,
140 double* u_grad_test_trace_ref,
141 double* normal_ref,
142 double* boundaryJac_ref,
143 //physics
144 int nElements_global,
145 int* u_l2g,
146 double* u_dof,
147 double alphaBDF,
148 double* q_vf,
149 double* q_divU,
150 double* q_vs,
151 double* q_vos,
152 double rho_s,
153 double* q_rho_f,
154 double rho_s_min,
155 double rho_f_min,
156 double* ebqe_vf,
157 double* ebqe_vs,
158 double* ebqe_vos,
159 double* ebqe_rho_f,
160 double* q_u,
161 double* q_grad_u,
162 double* ebqe_u,
163 double* ebqe_grad_u,
164 int offset_u,
165 int stride_u,
166 double* elementResidual_u,
167 int nExteriorElementBoundaries_global,
168 int* exteriorElementBoundariesArray,
169 int* elementBoundaryElementsArray,
170 int* elementBoundaryLocalElementBoundariesArray,
171 double* element_u,
172 int eN,
173 double compatibility_condition,
174 int INTEGRATE_BY_PARTS_DIV_U,
175 double* q_a)
176 {
177 for (int i=0;i<nDOF_test_element;i++)
178 {
179 elementResidual_u[i]=0.0;
180 }//i
181 //loop over quadrature points and compute integrands
182 for (int k=0;k<nQuadraturePoints_element;k++)
183 {
184 //compute indeces and declare local storage
185 int eN_k = eN*nQuadraturePoints_element+k,
186 eN_k_nSpace = eN_k*nSpace;
187 //eN_nDOF_trial_element = eN*nDOF_trial_element;
188 double u=0.0,grad_u[nSpace],
189 a=0.0,
190 f[nSpace],
191 jac[nSpace*nSpace],
192 jacDet,
193 jacInv[nSpace*nSpace],
194 u_grad_trial[nDOF_trial_element*nSpace],
195 u_test_dV[nDOF_trial_element],
196 u_grad_test_dV[nDOF_test_element*nSpace],
197 dV,x,y,z,
198 G[nSpace*nSpace],G_dd_G,tr_G;
199 //
200 //compute solution and gradients at quadrature points
201 //
202 ck.calculateMapping_element(eN,
203 k,
204 mesh_dof,
205 mesh_l2g,
206 mesh_trial_ref,
207 mesh_grad_trial_ref,
208 jac,
209 jacDet,
210 jacInv,
211 x,y,z);
212 //get the physical integration weight
213 dV = fabs(jacDet)*dV_ref[k];
214 ck.calculateG(jacInv,G,G_dd_G,tr_G);
215 //get the trial function gradients
216 ck.gradTrialFromRef(&u_grad_trial_ref[k*nDOF_trial_element*nSpace],jacInv,u_grad_trial);
217 //get the solution
218 ck.valFromElementDOF(element_u,&u_trial_ref[k*nDOF_trial_element],u);
219 //get the solution gradients
220 ck.gradFromElementDOF(element_u,u_grad_trial,grad_u);
221 //precalculate test function products with integration weights
222 for (int j=0;j<nDOF_trial_element;j++)
223 {
224 u_test_dV[j] = u_test_ref[k*nDOF_trial_element+j]*dV;
225 for (int I=0;I<nSpace;I++)
226 {
227 u_grad_test_dV[j*nSpace+I] = u_grad_trial[j*nSpace+I]*dV;//cek warning won't work for Petrov-Galerkin
228 }
229 }
230 //
231 //calculate pde coefficients at quadrature points
232 //
233 evaluateCoefficients(alphaBDF,
234 &q_vf[eN_k_nSpace],
235 &q_vs[eN_k_nSpace],
236 q_vos[eN_k],
237 rho_s_min,
238 rho_f_min,
239 f,
240 a);
241 //
242 //update element residual
243 //
244 for(int i=0;i<nDOF_test_element;i++)
245 {
246 //int eN_k_i=eN_k*nDOF_test_element+i;
247 //int eN_k_i_nSpace = eN_k_i*nSpace;
248 int i_nSpace=i*nSpace;
249 elementResidual_u[i] +=
250 (INTEGRATE_BY_PARTS_DIV_U == 1 ? ck.Advection_weak(f,&u_grad_test_dV[i_nSpace]) : q_divU[eN_k]*u_test_dV[i])
251 //+ compatibility_condition*u_test_dV[i] // mql: to make the system solvable if int(div(u))!=0
252 + ck.NumericalDiffusion(a,grad_u,&u_grad_test_dV[i_nSpace]);
253 }//i
254 //
255 //save momentum for time history and velocity for subgrid error
256 //save solution for other models
257 //
258 q_a[eN_k] = a;
259 q_u[eN_k] = u;
260
261 for (int I=0;I<nSpace;I++)
262 q_grad_u[eN_k_nSpace+I] = grad_u[I];
263 }
264 }
265
266 void calculateResidual(arguments_dict& args)
267 {
268 xt::pyarray<double>& mesh_trial_ref = args.array<double>("mesh_trial_ref");
269 xt::pyarray<double>& mesh_grad_trial_ref = args.array<double>("mesh_grad_trial_ref");
270 xt::pyarray<double>& mesh_dof = args.array<double>("mesh_dof");
271 xt::pyarray<int>& mesh_l2g = args.array<int>("mesh_l2g");
272 xt::pyarray<double>& dV_ref = args.array<double>("dV_ref");
273 xt::pyarray<double>& u_trial_ref = args.array<double>("u_trial_ref");
274 xt::pyarray<double>& u_grad_trial_ref = args.array<double>("u_grad_trial_ref");
275 xt::pyarray<double>& u_test_ref = args.array<double>("u_test_ref");
276 xt::pyarray<double>& u_grad_test_ref = args.array<double>("u_grad_test_ref");
277 xt::pyarray<double>& mesh_trial_trace_ref = args.array<double>("mesh_trial_trace_ref");
278 xt::pyarray<double>& mesh_grad_trial_trace_ref = args.array<double>("mesh_grad_trial_trace_ref");
279 xt::pyarray<double>& dS_ref = args.array<double>("dS_ref");
280 xt::pyarray<double>& u_trial_trace_ref = args.array<double>("u_trial_trace_ref");
281 xt::pyarray<double>& u_grad_trial_trace_ref = args.array<double>("u_grad_trial_trace_ref");
282 xt::pyarray<double>& u_test_trace_ref = args.array<double>("u_test_trace_ref");
283 xt::pyarray<double>& u_grad_test_trace_ref = args.array<double>("u_grad_test_trace_ref");
284 xt::pyarray<double>& normal_ref = args.array<double>("normal_ref");
285 xt::pyarray<double>& boundaryJac_ref = args.array<double>("boundaryJac_ref");
286 int nElements_global = args.scalar<int>("nElements_global");
287 xt::pyarray<int>& isDOFBoundary = args.array<int>("isDOFBoundary");
288 xt::pyarray<int>& isFluxBoundary = args.array<int>("isFluxBoundary");
289 xt::pyarray<int>& u_l2g = args.array<int>("u_l2g");
290 xt::pyarray<double>& u_dof = args.array<double>("u_dof");
291 double alphaBDF = args.scalar<double>("alphaBDF");
292 xt::pyarray<double>& q_vf = args.array<double>("q_vf");
293 xt::pyarray<double>& q_divU = args.array<double>("q_divU");
294 xt::pyarray<double>& q_vs = args.array<double>("q_vs");
295 xt::pyarray<double>& q_vos = args.array<double>("q_vos");
296 double rho_s = args.scalar<double>("rho_s");
297 xt::pyarray<double>& q_rho_f = args.array<double>("q_rho_f");
298 double rho_s_min = args.scalar<double>("rho_s_min");
299 double rho_f_min = args.scalar<double>("rho_f_min");
300 xt::pyarray<double>& ebqe_vf = args.array<double>("ebqe_vf");
301 xt::pyarray<double>& ebqe_vs = args.array<double>("ebqe_vs");
302 xt::pyarray<double>& ebqe_vos = args.array<double>("ebqe_vos");
303 xt::pyarray<double>& ebqe_rho_f = args.array<double>("ebqe_rho_f");
304 xt::pyarray<double>& q_u = args.array<double>("q_u");
305 xt::pyarray<double>& q_grad_u = args.array<double>("q_grad_u");
306 xt::pyarray<double>& ebqe_u = args.array<double>("ebqe_u");
307 xt::pyarray<double>& ebqe_grad_u = args.array<double>("ebqe_grad_u");
308 xt::pyarray<double>& ebqe_bc_u_ext = args.array<double>("ebqe_bc_u_ext");
309 xt::pyarray<double>& ebqe_adv_flux = args.array<double>("ebqe_adv_flux");
310 xt::pyarray<double>& ebqe_diff_flux = args.array<double>("ebqe_diff_flux");
311 xt::pyarray<double>& bc_adv_flux = args.array<double>("bc_adv_flux");
312 xt::pyarray<double>& bc_diff_flux = args.array<double>("bc_diff_flux");
313 int offset_u = args.scalar<int>("offset_u");
314 int stride_u = args.scalar<int>("stride_u");
315 xt::pyarray<double>& globalResidual = args.array<double>("globalResidual");
316 int nExteriorElementBoundaries_global = args.scalar<int>("nExteriorElementBoundaries_global");
317 xt::pyarray<int>& exteriorElementBoundariesArray = args.array<int>("exteriorElementBoundariesArray");
318 xt::pyarray<int>& elementBoundaryElementsArray = args.array<int>("elementBoundaryElementsArray");
319 xt::pyarray<int>& elementBoundaryLocalElementBoundariesArray = args.array<int>("elementBoundaryLocalElementBoundariesArray");
320 int INTEGRATE_BY_PARTS_DIV_U = args.scalar<int>("INTEGRATE_BY_PARTS_DIV_U");
321 xt::pyarray<double>& q_a = args.array<double>("q_a");
322 xt::pyarray<double>& ebqe_a = args.array<double>("ebqe_a");
323 double compatibility_condition=0.;
324 /* // COMPUTE COMPATIBILITY CONSTANT */
325 /* // mql: Modify the rhs (by adding a constant) so that the Poission system is solvable (assume diffusive flux = 0). */
326 /* // Note that this is equivalent to consider a (not known) diffusive flux != 0 s.t. the system is solvable. */
327 /* for(int eN=0;eN<nElements_global;eN++) */
328 /* { */
329 /* for (int k=0;k<nQuadraturePoints_element;k++) */
330 /* { */
331 /* int eN_k = eN*nQuadraturePoints_element+k; */
332 /* double */
333 /* jac[nSpace*nSpace], */
334 /* jacDet, */
335 /* jacInv[nSpace*nSpace], */
336 /* dV,x,y,z; */
337 /* ck.calculateMapping_element(eN, */
338 /* k, */
339 /* mesh_dof.data(), */
340 /* mesh_l2g.data(), */
341 /* mesh_trial_ref.data(), */
342 /* mesh_grad_trial_ref.data(), */
343 /* jac, */
344 /* jacDet, */
345 /* jacInv, */
346 /* x,y,z); */
347 /* //get the physical integration weight */
348 /* dV = fabs(jacDet)*dV_ref.data()[k]; */
349 /* compatibility_condition -= q_divU.data()[eN_k]*dV; */
350 /* } */
351 /* } */
352 //
353 //loop over elements to compute volume integrals and load them into element and global residual
354 //
355 //eN is the element index
356 //eN_k is the quadrature point index for a scalar
357 //eN_k_nSpace is the quadrature point index for a vector
358 //eN_i is the element test function index
359 //eN_j is the element trial function index
360 //eN_k_j is the quadrature point index for a trial function
361 //eN_k_i is the quadrature point index for a trial function
362 for(int eN=0;eN<nElements_global;eN++)
363 {
364 //declare local storage for element residual and initialize
365 double elementResidual_u[nDOF_test_element],element_u[nDOF_trial_element];
366 for (int i=0;i<nDOF_test_element;i++)
367 {
368 int eN_i=eN*nDOF_test_element+i;
369 element_u[i] = u_dof.data()[u_l2g.data()[eN_i]];
370 }//i
371 calculateElementResidual(mesh_trial_ref.data(),
372 mesh_grad_trial_ref.data(),
373 mesh_dof.data(),
374 mesh_l2g.data(),
375 dV_ref.data(),
376 u_trial_ref.data(),
377 u_grad_trial_ref.data(),
378 u_test_ref.data(),
379 u_grad_test_ref.data(),
380 mesh_trial_trace_ref.data(),
381 mesh_grad_trial_trace_ref.data(),
382 dS_ref.data(),
383 u_trial_trace_ref.data(),
384 u_grad_trial_trace_ref.data(),
385 u_test_trace_ref.data(),
386 u_grad_test_trace_ref.data(),
387 normal_ref.data(),
388 boundaryJac_ref.data(),
389 nElements_global,
390 u_l2g.data(),
391 u_dof.data(),
392 alphaBDF,
393 q_vf.data(),
394 q_divU.data(),
395 q_vs.data(),
396 q_vos.data(),
397 rho_s,
398 q_rho_f.data(),
399 rho_s_min,
400 rho_f_min,
401 ebqe_vf.data(),
402 ebqe_vs.data(),
403 ebqe_vos.data(),
404 ebqe_rho_f.data(),
405 q_u.data(),
406 q_grad_u.data(),
407 ebqe_u.data(),
408 ebqe_grad_u.data(),
409 offset_u,
410 stride_u,
411 elementResidual_u,
412 nExteriorElementBoundaries_global,
413 exteriorElementBoundariesArray.data(),
414 elementBoundaryElementsArray.data(),
415 elementBoundaryLocalElementBoundariesArray.data(),
416 element_u,
417 eN,
418 compatibility_condition,
419 INTEGRATE_BY_PARTS_DIV_U,
420 q_a.data());
421 //
422 //load element into global residual and save element residual
423 //
424 for(int i=0;i<nDOF_test_element;i++)
425 {
426 int eN_i=eN*nDOF_test_element+i;
427 globalResidual.data()[offset_u+stride_u*u_l2g.data()[eN_i]]+=elementResidual_u[i];
428 }//i
429 }//elements
430 //
431 //loop over exterior element boundaries to calculate levelset gradient
432 //
433 //ebNE is the Exterior element boundary INdex
434 //ebN is the element boundary INdex
435 //eN is the element index
436 for (int ebNE = 0; ebNE < nExteriorElementBoundaries_global; ebNE++)
437 {
438 int ebN = exteriorElementBoundariesArray.data()[ebNE],
439 eN = elementBoundaryElementsArray.data()[ebN*2+0],
440 ebN_local = elementBoundaryLocalElementBoundariesArray.data()[ebN*2+0];
441 //eN_nDOF_trial_element = eN*nDOF_trial_element;
442 double elementResidual_u[nDOF_test_element];
443 double element_u[nDOF_trial_element];
444 for (int i=0;i<nDOF_test_element;i++)
445 {
446 int eN_i=eN*nDOF_test_element+i;
447 element_u[i] = u_dof.data()[u_l2g.data()[eN_i]];
448 elementResidual_u[i] = 0.0;
449 }//i
450 for (int kb=0;kb<nQuadraturePoints_elementBoundary;kb++)
451 {
452 int ebNE_kb = ebNE*nQuadraturePoints_elementBoundary+kb,
453 ebNE_kb_nSpace = ebNE_kb*nSpace,
454 ebN_local_kb = ebN_local*nQuadraturePoints_elementBoundary+kb,
455 ebN_local_kb_nSpace = ebN_local_kb*nSpace;
456 double h_b=0.0,
457 penalty=0.0,
458 u_ext=0.0,
459 bc_u_ext=0.0,
460 adv_flux_ext=0.0,
461 diff_flux_ext=0.0,
462 a_ext=0.0,
463 f_ext[nSpace],
464 grad_u_ext[nSpace],
465 jac_ext[nSpace*nSpace],
466 jacDet_ext,
467 jacInv_ext[nSpace*nSpace],
468 boundaryJac[nSpace*(nSpace-1)],
469 metricTensor[(nSpace-1)*(nSpace-1)],
470 metricTensorDetSqrt,
471 dS,
472 u_test_dS[nDOF_test_element],
473 u_grad_trial_trace[nDOF_trial_element*nSpace],
474 u_grad_test_dS[nDOF_test_element*nSpace],
475 normal[nSpace],x_ext,y_ext,z_ext,
476 G[nSpace*nSpace],G_dd_G,tr_G;
477 //
478 //calculate the solution and gradients at quadrature points
479 //
480 ck.calculateMapping_elementBoundary(eN,
481 ebN_local,
482 kb,
483 ebN_local_kb,
484 mesh_dof.data(),
485 mesh_l2g.data(),
486 mesh_trial_trace_ref.data(),
487 mesh_grad_trial_trace_ref.data(),
488 boundaryJac_ref.data(),
489 jac_ext,
490 jacDet_ext,
491 jacInv_ext,
492 boundaryJac,
493 metricTensor,
494 metricTensorDetSqrt,
495 normal_ref.data(),
496 normal,
497 x_ext,y_ext,z_ext);
498 dS = metricTensorDetSqrt*dS_ref.data()[kb];
499 //get the metric tensor
500 //cek todo use symmetry
501 ck.calculateG(jacInv_ext,G,G_dd_G,tr_G);
502 ck.calculateGScale(G,normal,h_b);
503 penalty = 10.0/h_b;
504 //compute shape and solution information
505 //shape
506 ck.gradTrialFromRef(&u_grad_trial_trace_ref.data()[ebN_local_kb_nSpace*nDOF_trial_element],jacInv_ext,u_grad_trial_trace);
507 //solution and gradients
508 ck.valFromElementDOF(element_u,&u_trial_trace_ref.data()[ebN_local_kb*nDOF_test_element],u_ext);
509 ck.gradFromElementDOF(element_u,u_grad_trial_trace,grad_u_ext);
510 //precalculate test function products with integration weights
511 for (int j=0;j<nDOF_trial_element;j++)
512 {
513 u_test_dS[j] = u_test_trace_ref.data()[ebN_local_kb*nDOF_test_element+j]*dS;
514 for (int I=0;I<nSpace;I++)
515 u_grad_test_dS[j*nSpace+I] = u_grad_trial_trace[j*nSpace+I]*dS;//cek hack, using trial
516 }
517 //
518 //load the boundary values
519 //
520 bc_u_ext = isDOFBoundary.data()[ebNE_kb]*ebqe_bc_u_ext.data()[ebNE_kb]+(1-isDOFBoundary.data()[ebNE_kb])*u_ext;
521 //
522 //calculate the pde coefficients using the solution and the boundary values for the solution
523 //
524
525 evaluateCoefficients(alphaBDF,
526 &ebqe_vf.data()[ebNE_kb_nSpace],
527 &ebqe_vs.data()[ebNE_kb_nSpace],
528 q_vos.data()[ebNE_kb],
529 rho_s_min,
530 rho_f_min,
531 f_ext,
532 a_ext);
533 ebqe_u.data()[ebNE_kb] = u_ext;
534 for (int I=0;I<nSpace;I++)
535 ebqe_grad_u.data()[ebNE_kb_nSpace+I] = grad_u_ext[I];
536 //
537 //calculate the numerical fluxes
538 //
539 exteriorNumericalAdvectiveFlux(isFluxBoundary.data()[ebNE_kb],
540 bc_adv_flux.data()[ebNE_kb],
541 normal,
542 f_ext,
543 adv_flux_ext); //=f.normal = [(1-vos)*vf + vos*vs].normal
544 exteriorNumericalDiffusiveFlux(isDOFBoundary.data()[ebNE_kb],
545 isFluxBoundary.data()[ebNE_kb],
546 normal,
547 a_ext,
548 grad_u_ext,
549 u_ext,
550 bc_u_ext,
551 bc_diff_flux.data()[ebNE_kb],
552 penalty,
553 diff_flux_ext);
554 if (isDOFBoundary.data()[ebNE_kb] != 1)
555 diff_flux_ext = 0.0; // mql: don't consider diffusive flux unless Dirichlet BC
556 //if(isFluxBoundary.data()[ebNE_kb] == 1)
557 //{
558 // adv_flux_ext = 0.0;
559 // diff_flux_ext = 0.0;
560 //}
561 ebqe_a.data()[ebNE_kb] = a_ext;
562 ebqe_adv_flux.data()[ebNE_kb] = adv_flux_ext;
563 ebqe_diff_flux.data()[ebNE_kb] = diff_flux_ext;
564 //
565 //update residuals
566 //
567 for (int i=0;i<nDOF_test_element;i++)
568 {
569 elementResidual_u[i] +=
570 (INTEGRATE_BY_PARTS_DIV_U == 1 ? ck.ExteriorElementBoundaryFlux(adv_flux_ext,u_test_dS[i]) : 0.)
571 + ck.ExteriorElementBoundaryFlux(diff_flux_ext,u_test_dS[i]) // mql: just != 0 if Dirichlet BC
572 + ck.ExteriorElementBoundaryScalarDiffusionAdjoint(isDOFBoundary.data()[ebNE_kb],
573 0, // mql: if Dirichlet BCs, then the flux BCs don't matter
574 //isFluxBoundary.data()[ebNE_kb],
575 1.0,
576 u_ext,
577 bc_u_ext,
578 normal,
579 a_ext,
580 &u_grad_test_dS[i*nSpace]);
581 }//i
582 }//kb
583 //
584 //update the element and global residual storage
585 //
586 for (int i=0;i<nDOF_test_element;i++)
587 {
588 int eN_i = eN*nDOF_test_element+i;
589 globalResidual.data()[offset_u+stride_u*u_l2g.data()[eN_i]] += elementResidual_u[i];
590 }//i
591 }//ebNE
592
593 }
594
595 inline void calculateElementJacobian(//element
596 double* mesh_trial_ref,
597 double* mesh_grad_trial_ref,
598 double* mesh_dof,
599 int* mesh_l2g,
600 double* dV_ref,
601 double* u_trial_ref,
602 double* u_grad_trial_ref,
603 double* u_test_ref,
604 double* u_grad_test_ref,
605 //element boundary
606 double* mesh_trial_trace_ref,
607 double* mesh_grad_trial_trace_ref,
608 double* dS_ref,
609 double* u_trial_trace_ref,
610 double* u_grad_trial_trace_ref,
611 double* u_test_trace_ref,
612 double* u_grad_test_trace_ref,
613 double* normal_ref,
614 double* boundaryJac_ref,
615 //physics
616 int nElements_global,
617 int* u_l2g,
618 double* u_dof,
619 double alphaBDF,
620 double* q_vf,
621 double* q_vs,
622 double* q_vos,
623 double rho_s,
624 double* q_rho_f,
625 double rho_s_min,
626 double rho_f_min,
627 double* elementJacobian_u_u,
628 double* element_u,
629 int eN)
630 {
631 for (int i=0;i<nDOF_test_element;i++)
632 for (int j=0;j<nDOF_trial_element;j++)
633 {
634 elementJacobian_u_u[i*nDOF_trial_element+j]=0.0;
635 }
636 for (int k=0;k<nQuadraturePoints_element;k++)
637 {
638 int eN_k = eN*nQuadraturePoints_element+k, //index to a scalar at a quadrature point
639 eN_k_nSpace = eN_k*nSpace;
640 //eN_nDOF_trial_element = eN*nDOF_trial_element; //index to a vector at a quadrature point
641
642 //declare local storage
643 double u=0.0,
644 grad_u[nSpace],
645 f[nSpace],
646 a=0.0,
647 jac[nSpace*nSpace],
648 jacDet,
649 jacInv[nSpace*nSpace],
650 u_grad_trial[nDOF_trial_element*nSpace],
651 dV,
652 u_grad_test_dV[nDOF_test_element*nSpace],
653 x,y,z,
654 G[nSpace*nSpace],G_dd_G,tr_G;
655 //
656 //calculate solution and gradients at quadrature points
657 //
658 ck.calculateMapping_element(eN,
659 k,
660 mesh_dof,
661 mesh_l2g,
662 mesh_trial_ref,
663 mesh_grad_trial_ref,
664 jac,
665 jacDet,
666 jacInv,
667 x,y,z);
668 //get the physical integration weight
669 dV = fabs(jacDet)*dV_ref[k];
670 ck.calculateG(jacInv,G,G_dd_G,tr_G);
671 //get the trial function gradients
672 ck.gradTrialFromRef(&u_grad_trial_ref[k*nDOF_trial_element*nSpace],jacInv,u_grad_trial);
673 //get the solution
674 ck.valFromElementDOF(element_u,&u_trial_ref[k*nDOF_trial_element],u);
675 //get the solution gradients
676 ck.gradFromElementDOF(element_u,u_grad_trial,grad_u);
677 //precalculate test function products with integration weights
678 for (int j=0;j<nDOF_trial_element;j++)
679 {
680 for (int I=0;I<nSpace;I++)
681 {
682 u_grad_test_dV[j*nSpace+I] = u_grad_trial[j*nSpace+I]*dV;//cek warning won't work for Petrov-Galerkin
683 }
684 }
685 //
686 //calculate pde coefficients and derivatives at quadrature points
687 //
688 evaluateCoefficients(alphaBDF,
689 &q_vf[eN_k_nSpace],
690 &q_vs[eN_k_nSpace],
691 q_vos[eN_k],
692 rho_s_min,
693 rho_f_min,
694 f,
695 a);
696 for(int i=0;i<nDOF_test_element;i++)
697 {
698 //int eN_k_i=eN_k*nDOF_test_element+i;
699 //int eN_k_i_nSpace=eN_k_i*nSpace;
700 int i_nSpace=i*nSpace;
701 for(int j=0;j<nDOF_trial_element;j++)
702 {
703 //int eN_k_j=eN_k*nDOF_trial_element+j;
704 //int eN_k_j_nSpace = eN_k_j*nSpace;
705 int j_nSpace = j*nSpace;
706 elementJacobian_u_u[i*nDOF_trial_element+j] +=
707 ck.NumericalDiffusionJacobian(a,&u_grad_trial[j_nSpace],&u_grad_test_dV[i_nSpace]);
708 }//j
709 }//i
710 }//k
711 }
712
713 void calculateJacobian(arguments_dict& args)
714 {
715 xt::pyarray<double>& mesh_trial_ref = args.array<double>("mesh_trial_ref");
716 xt::pyarray<double>& mesh_grad_trial_ref = args.array<double>("mesh_grad_trial_ref");
717 xt::pyarray<double>& mesh_dof = args.array<double>("mesh_dof");
718 xt::pyarray<int>& mesh_l2g = args.array<int>("mesh_l2g");
719 xt::pyarray<double>& dV_ref = args.array<double>("dV_ref");
720 xt::pyarray<double>& u_trial_ref = args.array<double>("u_trial_ref");
721 xt::pyarray<double>& u_grad_trial_ref = args.array<double>("u_grad_trial_ref");
722 xt::pyarray<double>& u_test_ref = args.array<double>("u_test_ref");
723 xt::pyarray<double>& u_grad_test_ref = args.array<double>("u_grad_test_ref");
724 xt::pyarray<double>& mesh_trial_trace_ref = args.array<double>("mesh_trial_trace_ref");
725 xt::pyarray<double>& mesh_grad_trial_trace_ref = args.array<double>("mesh_grad_trial_trace_ref");
726 xt::pyarray<double>& dS_ref = args.array<double>("dS_ref");
727 xt::pyarray<double>& u_trial_trace_ref = args.array<double>("u_trial_trace_ref");
728 xt::pyarray<double>& u_grad_trial_trace_ref = args.array<double>("u_grad_trial_trace_ref");
729 xt::pyarray<double>& u_test_trace_ref = args.array<double>("u_test_trace_ref");
730 xt::pyarray<double>& u_grad_test_trace_ref = args.array<double>("u_grad_test_trace_ref");
731 xt::pyarray<double>& normal_ref = args.array<double>("normal_ref");
732 xt::pyarray<double>& boundaryJac_ref = args.array<double>("boundaryJac_ref");
733 int nElements_global = args.scalar<int>("nElements_global");
734 xt::pyarray<int>& isDOFBoundary = args.array<int>("isDOFBoundary");
735 xt::pyarray<int>& isFluxBoundary = args.array<int>("isFluxBoundary");
736 xt::pyarray<int>& u_l2g = args.array<int>("u_l2g");
737 xt::pyarray<double>& u_dof = args.array<double>("u_dof");
738 double alphaBDF = args.scalar<double>("alphaBDF");
739 xt::pyarray<double>& q_vf = args.array<double>("q_vf");
740 xt::pyarray<double>& q_vs = args.array<double>("q_vs");
741 xt::pyarray<double>& q_vos = args.array<double>("q_vos");
742 double rho_s = args.scalar<double>("rho_s");
743 xt::pyarray<double>& q_rho_f = args.array<double>("q_rho_f");
744 double rho_s_min = args.scalar<double>("rho_s_min");
745 double rho_f_min = args.scalar<double>("rho_f_min");
746 xt::pyarray<double>& ebqe_vf = args.array<double>("ebqe_vf");
747 xt::pyarray<double>& ebqe_vs = args.array<double>("ebqe_vs");
748 xt::pyarray<double>& ebqe_vos = args.array<double>("ebqe_vos");
749 xt::pyarray<double>& ebqe_rho_f = args.array<double>("ebqe_rho_f");
750 xt::pyarray<int>& csrRowIndeces_u_u = args.array<int>("csrRowIndeces_u_u");
751 xt::pyarray<int>& csrColumnOffsets_u_u = args.array<int>("csrColumnOffsets_u_u");
752 xt::pyarray<double>& globalJacobian = args.array<double>("globalJacobian");
753 int nExteriorElementBoundaries_global = args.scalar<int>("nExteriorElementBoundaries_global");
754 xt::pyarray<int>& exteriorElementBoundariesArray = args.array<int>("exteriorElementBoundariesArray");
755 xt::pyarray<int>& elementBoundaryElementsArray = args.array<int>("elementBoundaryElementsArray");
756 xt::pyarray<int>& elementBoundaryLocalElementBoundariesArray = args.array<int>("elementBoundaryLocalElementBoundariesArray");
757 xt::pyarray<int>& csrColumnOffsets_eb_u_u = args.array<int>("csrColumnOffsets_eb_u_u");
758 //
759 //loop over elements to compute volume integrals and load them into the element Jacobians and global Jacobian
760 //
761 for(int eN=0;eN<nElements_global;eN++)
762 {
763 double elementJacobian_u_u[nDOF_test_element*nDOF_trial_element],element_u[nDOF_trial_element];
764 for (int j=0;j<nDOF_trial_element;j++)
765 {
766 int eN_j = eN*nDOF_trial_element+j;
767 element_u[j] = u_dof.data()[u_l2g.data()[eN_j]];
768 }
769 calculateElementJacobian(mesh_trial_ref.data(),
770 mesh_grad_trial_ref.data(),
771 mesh_dof.data(),
772 mesh_l2g.data(),
773 dV_ref.data(),
774 u_trial_ref.data(),
775 u_grad_trial_ref.data(),
776 u_test_ref.data(),
777 u_grad_test_ref.data(),
778 mesh_trial_trace_ref.data(),
779 mesh_grad_trial_trace_ref.data(),
780 dS_ref.data(),
781 u_trial_trace_ref.data(),
782 u_grad_trial_trace_ref.data(),
783 u_test_trace_ref.data(),
784 u_grad_test_trace_ref.data(),
785 normal_ref.data(),
786 boundaryJac_ref.data(),
787 nElements_global,
788 u_l2g.data(),
789 u_dof.data(),
790 alphaBDF,
791 q_vf.data(),
792 q_vs.data(),
793 q_vos.data(),
794 rho_s,
795 q_rho_f.data(),
796 rho_s_min,
797 rho_f_min,
798 elementJacobian_u_u,
799 element_u,
800 eN);
801 //
802 //load into element Jacobian into global Jacobian
803 //
804 for (int i=0;i<nDOF_test_element;i++)
805 {
806 int eN_i = eN*nDOF_test_element+i;
807 for (int j=0;j<nDOF_trial_element;j++)
808 {
809 int eN_i_j = eN_i*nDOF_trial_element+j;
810 globalJacobian.data()[csrRowIndeces_u_u.data()[eN_i] + csrColumnOffsets_u_u.data()[eN_i_j]] += elementJacobian_u_u[i*nDOF_trial_element+j];
811 }//j
812 }//i
813 }//elements
814 //
815 //loop over exterior element boundaries to compute the surface integrals and load them into the global Jacobian
816 //
817 for (int ebNE = 0; ebNE < nExteriorElementBoundaries_global; ebNE++)
818 {
819 int ebN = exteriorElementBoundariesArray.data()[ebNE];
820 int eN = elementBoundaryElementsArray.data()[ebN*2+0],
821 ebN_local = elementBoundaryLocalElementBoundariesArray.data()[ebN*2+0],
822 eN_nDOF_trial_element = eN*nDOF_trial_element;
823 for (int kb=0;kb<nQuadraturePoints_elementBoundary;kb++)
824 {
825 int ebNE_kb = ebNE*nQuadraturePoints_elementBoundary+kb,
826 ebNE_kb_nSpace = ebNE_kb*nSpace,
827 ebN_local_kb = ebN_local*nQuadraturePoints_elementBoundary+kb,
828 ebN_local_kb_nSpace = ebN_local_kb*nSpace;
829 double h_b=0.0,
830 u_ext=0.0,
831 grad_u_ext[nSpace],
832 a_ext=0.0,
833 f_ext[nSpace],
834 jac_ext[nSpace*nSpace],
835 jacDet_ext,
836 jacInv_ext[nSpace*nSpace],
837 boundaryJac[nSpace*(nSpace-1)],
838 metricTensor[(nSpace-1)*(nSpace-1)],
839 metricTensorDetSqrt,
840 dS,
841 u_test_dS[nDOF_test_element],
842 u_grad_trial_trace[nDOF_trial_element*nSpace],
843 u_grad_test_dS[nDOF_test_element*nSpace],
844 normal[nSpace],x_ext,y_ext,z_ext,
845 penalty=0.0,
846 //
847 G[nSpace*nSpace],G_dd_G,tr_G;
848 //
849 //calculate the solution and gradients at quadrature points
850 //
851 u_ext=0.0;
852 for (int I=0;I<nSpace;I++)
853 {
854 grad_u_ext[I] = 0.0;
855 }
856 ck.calculateMapping_elementBoundary(eN,
857 ebN_local,
858 kb,
859 ebN_local_kb,
860 mesh_dof.data(),
861 mesh_l2g.data(),
862 mesh_trial_trace_ref.data(),
863 mesh_grad_trial_trace_ref.data(),
864 boundaryJac_ref.data(),
865 jac_ext,
866 jacDet_ext,
867 jacInv_ext,
868 boundaryJac,
869 metricTensor,
870 metricTensorDetSqrt,
871 normal_ref.data(),
872 normal,
873 x_ext,y_ext,z_ext);
874 dS = metricTensorDetSqrt*dS_ref.data()[kb];
875 ck.calculateG(jacInv_ext,G,G_dd_G,tr_G);
876 ck.calculateGScale(G,normal,h_b);
877 penalty=10.0/h_b;
878 //compute shape and solution information
879 //shape
880 ck.gradTrialFromRef(&u_grad_trial_trace_ref.data()[ebN_local_kb_nSpace*nDOF_trial_element],jacInv_ext,u_grad_trial_trace);
881 //solution and gradients
882 ck.valFromDOF(u_dof.data(),&u_l2g.data()[eN_nDOF_trial_element],&u_trial_trace_ref.data()[ebN_local_kb*nDOF_test_element],u_ext);
883 ck.gradFromDOF(u_dof.data(),&u_l2g.data()[eN_nDOF_trial_element],u_grad_trial_trace,grad_u_ext);
884 //precalculate test function products with integration weights
885 for (int j=0;j<nDOF_trial_element;j++)
886 {
887 u_test_dS[j] = u_test_trace_ref.data()[ebN_local_kb*nDOF_test_element+j]*dS;
888 for (int I=0;I<nSpace;I++)
889 u_grad_test_dS[j*nSpace+I] = u_grad_trial_trace[j*nSpace+I]*dS;//cek hack, using trial
890 }
891 //
892 //calculate the internal and external trace of the pde coefficients
893 //
894 evaluateCoefficients(alphaBDF,
895 &ebqe_vf.data()[ebNE_kb_nSpace],
896 &ebqe_vs.data()[ebNE_kb_nSpace],
897 q_vos.data()[ebNE_kb],
898 rho_s_min,
899 rho_f_min,
900 f_ext,
901 a_ext);
902 //
903 //update the global Jacobian from the flux Jacobian
904 //
905 for (int i=0;i<nDOF_test_element;i++)
906 {
907 int eN_i = eN*nDOF_test_element+i;
908 //int ebNE_kb_i = ebNE_kb*nDOF_test_element+i;
909 for (int j=0;j<nDOF_trial_element;j++)
910 {
911 int ebN_i_j = ebN*4*nDOF_test_X_trial_element + i*nDOF_trial_element + j,
912 ebN_local_kb_j=ebN_local_kb*nDOF_trial_element+j,
913 j_nSpace = j*nSpace;
914
915 globalJacobian.data()[csrRowIndeces_u_u.data()[eN_i] + csrColumnOffsets_eb_u_u.data()[ebN_i_j]] +=
916 ExteriorNumericalDiffusiveFluxJacobian(isDOFBoundary.data()[ebNE_kb],
917 isFluxBoundary.data()[ebNE_kb],
918 normal,
919 a_ext,
920 u_trial_trace_ref.data()[ebN_local_kb_j],
921 &u_grad_trial_trace[j_nSpace],
922 penalty)*u_test_dS[i]
923 +
924 ck.ExteriorElementBoundaryScalarDiffusionAdjointJacobian
925 (isDOFBoundary.data()[ebNE_kb],
926 isFluxBoundary.data()[ebNE_kb],
927 1.0,
928 u_trial_trace_ref.data()[ebN_local_kb_j],
929 normal,
930 a_ext,
931 &u_grad_test_dS[i*nSpace]);
932 }//j
933 }//i
934 }//kb
935 }//ebNE
936 }//computeJacobian
937 };//cppPresInc
938
939 inline cppPresInc_base* newPresInc(int nSpaceIn,
940 int nQuadraturePoints_elementIn,
941 int nDOF_mesh_trial_elementIn,
942 int nDOF_trial_elementIn,
943 int nDOF_test_elementIn,
944 int nQuadraturePoints_elementBoundaryIn,
945 int CompKernelFlag)
946 {
947 if (nSpaceIn == 2)
948 return proteus::chooseAndAllocateDiscretization2D<cppPresInc_base,cppPresInc,CompKernel>(nSpaceIn,
949 nQuadraturePoints_elementIn,
950 nDOF_mesh_trial_elementIn,
951 nDOF_trial_elementIn,
952 nDOF_test_elementIn,
953 nQuadraturePoints_elementBoundaryIn,
954 CompKernelFlag);
955 else
956 return proteus::chooseAndAllocateDiscretization<cppPresInc_base,cppPresInc,CompKernel>(nSpaceIn,
957 nQuadraturePoints_elementIn,
958 nDOF_mesh_trial_elementIn,
959 nDOF_trial_elementIn,
960 nDOF_test_elementIn,
961 nQuadraturePoints_elementBoundaryIn,
962 CompKernelFlag);
963 }
964}//proteus
965#endif
n
Int n
Definition Headers.h:28
u
Double u
Definition Headers.h:89
z
Double * z
Definition Headers.h:49
v
Double v
Definition Headers.h:95
proteus::cppPresInc_base
Definition PresInc.h:15
proteus::cppPresInc_base::~cppPresInc_base
virtual ~cppPresInc_base()
Definition PresInc.h:17
proteus::cppPresInc_base::calculateJacobian
virtual void calculateJacobian(arguments_dict &args)=0
proteus::cppPresInc_base::calculateResidual
virtual void calculateResidual(arguments_dict &args)=0
proteus::cppPresInc::ck
CompKernelType ck
Definition PresInc.h:33
proteus::cppPresInc::calculateElementJacobian
void calculateElementJacobian(double *mesh_trial_ref, double *mesh_grad_trial_ref, double *mesh_dof, int *mesh_l2g, double *dV_ref, double *u_trial_ref, double *u_grad_trial_ref, double *u_test_ref, double *u_grad_test_ref, double *mesh_trial_trace_ref, double *mesh_grad_trial_trace_ref, double *dS_ref, double *u_trial_trace_ref, double *u_grad_trial_trace_ref, double *u_test_trace_ref, double *u_grad_test_trace_ref, double *normal_ref, double *boundaryJac_ref, int nElements_global, int *u_l2g, double *u_dof, double alphaBDF, double *q_vf, double *q_vs, double *q_vos, double rho_s, double *q_rho_f, double rho_s_min, double rho_f_min, double *elementJacobian_u_u, double *element_u, int eN)
Definition PresInc.h:595
proteus::cppPresInc::calculateJacobian
void calculateJacobian(arguments_dict &args)
Definition PresInc.h:713
proteus::cppPresInc::exteriorNumericalDiffusiveFlux
void exteriorNumericalDiffusiveFlux(const int &isDOFBoundary, const int &isFluxBoundary, const double n[nSpace], const double &a, const double grad_potential[nSpace], const double &u, const double &bc_u, const double &bc_flux, const double &penalty, double &flux)
Definition PresInc.h:75
proteus::cppPresInc::ExteriorNumericalDiffusiveFluxJacobian
double ExteriorNumericalDiffusiveFluxJacobian(const int &isDOFBoundary, const int &isFluxBoundary, const double n[nSpace], const double &a, const double &v, const double grad_v[nSpace], const double &penalty)
Definition PresInc.h:105
proteus::cppPresInc::evaluateCoefficients
void evaluateCoefficients(const double &alphaBDF, const double vf[nSpace], const double vs[nSpace], const double &vos, const double &rhos_min, const double &rhof_min, double f[nSpace], double &a)
Definition PresInc.h:39
proteus::cppPresInc::exteriorNumericalAdvectiveFlux
void exteriorNumericalAdvectiveFlux(const int &isFluxBoundary, const double &bc_flux, const double n[nSpace], const double f[nSpace], double &flux)
Definition PresInc.h:58
proteus::cppPresInc::calculateElementResidual
void calculateElementResidual(double *mesh_trial_ref, double *mesh_grad_trial_ref, double *mesh_dof, int *mesh_l2g, double *dV_ref, double *u_trial_ref, double *u_grad_trial_ref, double *u_test_ref, double *u_grad_test_ref, double *mesh_trial_trace_ref, double *mesh_grad_trial_trace_ref, double *dS_ref, double *u_trial_trace_ref, double *u_grad_trial_trace_ref, double *u_test_trace_ref, double *u_grad_test_trace_ref, double *normal_ref, double *boundaryJac_ref, int nElements_global, int *u_l2g, double *u_dof, double alphaBDF, double *q_vf, double *q_divU, double *q_vs, double *q_vos, double rho_s, double *q_rho_f, double rho_s_min, double rho_f_min, double *ebqe_vf, double *ebqe_vs, double *ebqe_vos, double *ebqe_rho_f, double *q_u, double *q_grad_u, double *ebqe_u, double *ebqe_grad_u, int offset_u, int stride_u, double *elementResidual_u, int nExteriorElementBoundaries_global, int *exteriorElementBoundariesArray, int *elementBoundaryElementsArray, int *elementBoundaryLocalElementBoundariesArray, double *element_u, int eN, double compatibility_condition, int INTEGRATE_BY_PARTS_DIV_U, double *q_a)
Definition PresInc.h:123
proteus::cppPresInc::calculateResidual
void calculateResidual(arguments_dict &args)
Definition PresInc.h:266
proteus::cppPresInc::nDOF_test_X_trial_element
const int nDOF_test_X_trial_element
Definition PresInc.h:32
proteus
Definition ADR.h:19
proteus::newPresInc
cppPresInc_base * newPresInc(int nSpaceIn, int nQuadraturePoints_elementIn, int nDOF_mesh_trial_elementIn, int nDOF_trial_elementIn, int nDOF_test_elementIn, int nQuadraturePoints_elementBoundaryIn, int CompKernelFlag)
Definition PresInc.h:939
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::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::f
double f(const double &g, const double &h, const double &hZ)
Definition SW2DCV.h:58
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