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