proteus 1.9.0
C/C++/Fortran libraries
Loading...
Searching...
No Matches
MCorr.h
Go to the documentation of this file.
1#ifndef MCorr_H
2#define MCorr_H
3#include <cmath>
4#include <iostream>
5#include <set>
6#include <map>
7#include <valarray>
8#include "CompKernel.h"
9#include "ModelFactory.h"
10#include "equivalent_polynomials.h"
11#include PROTEUS_LAPACK_H
12#include "ArgumentsDict.h"
13#include "xtensor-python/pyarray.hpp"
14
15namespace py = pybind11;
16
17namespace proteus
18{
19
20 template<int nSpace, int nP, int nQ, int nEBQ>
21 // The trailing flag is the IFEM gate. It has always been false here -- by
22 // default rather than by statement -- and must stay false: with it set, an
23 // element whose interface passes through an edge or corner node takes the
24 // IFEM branch in Simplex::set_quad, which forces D to 0 and H/ImH to a hard
25 // 0/1 instead of the moment fit. That deletes the interface measure this
26 // model integrates over. Written out so the choice is visible at the call
27 // site and cannot change underneath us if the template default changes.
28 using GeneralizedFunctions = equivalent_polynomials::GeneralizedFunctions_mix<nSpace, nP, nP, nQ, nEBQ, false>;
29 //using GeneralizedFunctions = equivalent_polynomials::Regularized<nSpace, nP, nQ>;
30 //using GeneralizedFunctions = equivalent_polynomials::EquivalentPolynomials<nSpace, nP, nQ>;
31
32 class MCorr_base
33 {
34 public:
35 std::valarray<double> Rpos, Rneg;
36 std::valarray<double> FluxCorrectionMatrix;
37 virtual ~MCorr_base(){}
38 virtual void calculateResidual(arguments_dict& args, bool useExact)=0;
39 virtual void calculateJacobian(arguments_dict& args, bool useExact)=0;
40 virtual void elementSolve(arguments_dict& args)=0;
41 virtual void elementConstantSolve(arguments_dict& args)=0;
42 virtual std::tuple<double, double> globalConstantRJ(arguments_dict& args)=0;
43 virtual double calculateMass(arguments_dict& args, bool useExact)=0;
44 virtual void setMassQuadrature(arguments_dict& args, bool useExact)=0;
45 virtual void FCTStep(arguments_dict& args)=0;
46 virtual void calculateMassMatrix(arguments_dict& args)=0;
47 virtual void setMassQuadratureEdgeBasedStabilizationMethods(arguments_dict& args,
48 bool useExact)=0;
49 };
50
51 template<class CompKernelType,
52 int nSpace,
53 int nQuadraturePoints_element,
54 int nDOF_mesh_trial_element,
55 int nDOF_trial_element,
56 int nDOF_test_element,
57 int nQuadraturePoints_elementBoundary>
58 class MCorr : public MCorr_base
59 {
60 public:
61 std::set<int> cutfem_boundaries;
62 std::map<int, int> cutfem_local_boundaries;
63 const int nDOF_test_X_trial_element;
64 CompKernelType ck;
65 GeneralizedFunctions<nSpace,2,nQuadraturePoints_element,nQuadraturePoints_elementBoundary> gf;
66 GeneralizedFunctions<nSpace,2,nDOF_trial_element,nQuadraturePoints_elementBoundary> gf_nodes;
67 GeneralizedFunctions<nSpace,4,nQuadraturePoints_element,nQuadraturePoints_elementBoundary> gf_s;
68 MCorr():
69 nDOF_test_X_trial_element(nDOF_test_element*nDOF_trial_element),
70 ck()
71 {}
72
73 inline
74 void evaluateCoefficients(const double& epsHeaviside,
75 const double& epsDirac,
76 const double& phi,
77 const double& H,
78 const double& u,
79 const double& porosity,
80 double& r,
81 double& dr)
82 {
83 r = porosity*(gf.H(epsHeaviside,phi+u) - H);
84 dr = porosity*gf.D(epsDirac,phi+u);
85 }
86
87 inline void calculateElementResidual(//element
88 double* mesh_trial_ref,
89 double* mesh_grad_trial_ref,
90 double* mesh_dof,
91 int* mesh_l2g,
92 double* dV_ref,
93 double* u_trial_ref,
94 double* u_grad_trial_ref,
95 double* u_test_ref,
96 double* u_grad_test_ref,
97 //element boundary
98 double* mesh_trial_trace_ref,
99 double* mesh_grad_trial_trace_ref,
100 double* dS_ref,
101 double* u_trial_trace_ref,
102 double* u_grad_trial_trace_ref,
103 double* u_test_trace_ref,
104 double* u_grad_test_trace_ref,
105 double* normal_ref,
106 double* boundaryJac_ref,
107 //physics
108 int nElements_global,
109 double useMetrics,
110 double epsFactHeaviside,
111 double epsFactDirac,
112 double epsFactDiffusion,
113 int* u_l2g,
114 int* r_l2g,
115 double* elementDiameter,
116 double* nodeDiametersArray,
117 double* u_dof,
118 double* q_phi,
119 double* q_normal_phi,
120 double* ebqe_phi,
121 double* ebqe_normal_phi,
122 double* q_H,
123 double* q_u,
124 double* q_n,
125 double* ebqe_u,
126 double* ebqe_n,
127 double* q_r,
128 double* q_porosity,
129 int offset_u, int stride_u,
130 double* elementResidual_u,
131 int nExteriorElementBoundaries_global,
132 int* exteriorElementBoundariesArray,
133 int* elementBoundaryElementsArray,
134 int* elementBoundaryLocalElementBoundariesArray,
135 double* element_u,
136 int eN,
137 bool element_active,
138 double* isActiveR,
139 double* isActiveDOF,
140 const double* phi_solid)
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 gf.set_quad(k);
165 gf_s.set_quad(k);
166 const double H_s = gf_s.H(0.0,phi_solid[eN_k]);
167 //
168 //compute solution and gradients at quadrature points
169 //
170 ck.calculateMapping_element(eN,
171 k,
172 mesh_dof,
173 mesh_l2g,
174 mesh_trial_ref,
175 mesh_grad_trial_ref,
176 jac,
177 jacDet,
178 jacInv,
179 x,y,z);
180 ck.calculateH_element(eN,
181 k,
182 nodeDiametersArray,
183 mesh_l2g,
184 mesh_trial_ref,
185 h_phi);
186 //get the physical integration weight
187 dV = fabs(jacDet)*dV_ref[k];
188 ck.calculateG(jacInv,G,G_dd_G,tr_G);
189
190 /* double dir[nSpace]; */
191 /* double norm = 1.0e-8; */
192 /* for (int I=0;I<nSpace;I++) */
193 /* norm += q_normal_phi[eN_k_nSpace+I]*q_normal_phi[eN_k_nSpace+I]; */
194 /* norm = sqrt(norm); */
195 /* for (int I=0;I<nSpace;I++) */
196 /* dir[I] = q_normal_phi[eN_k_nSpace+I]/norm; */
197 /* ck.calculateGScale(G,dir,h_phi); */
198
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 u_test_dV[j] = u_test_ref[k*nDOF_trial_element+j]*dV;
209 for (int I=0;I<nSpace;I++)
210 {
211 u_grad_test_dV[j*nSpace+I] = u_grad_trial[j*nSpace+I]*dV;//cek warning won't work for Petrov-Galerkin
212 }
213 }
214
215
216
217 //
218 //calculate pde coefficients at quadrature points
219 //
220 epsHeaviside = epsFactHeaviside*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter[eN]);
221 epsDirac = epsFactDirac* (useMetrics*h_phi+(1.0-useMetrics)*elementDiameter[eN]);
222 epsDiffusion = epsFactDiffusion*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter[eN]);
223 // *(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter[eN]);
224 evaluateCoefficients(epsHeaviside,
225 epsDirac,
226 q_phi[eN_k],
227 q_H[eN_k],
228 u,
229 q_porosity[eN_k],
230 r,
231 dr);
232 //
233 //update element residual
234 //
235 for(int i=0;i<nDOF_test_element;i++)
236 {
237 int eN_i = eN*nDOF_test_element+i;
238 //int eN_k_i=eN_k*nDOF_test_element+i;
239 //int eN_k_i_nSpace = eN_k_i*nSpace;
240 int i_nSpace=i*nSpace;
241
242 elementResidual_u[i] += H_s*(ck.Reaction_weak(r,u_test_dV[i]) +
243 ck.NumericalDiffusion(epsDiffusion,grad_u,&u_grad_test_dV[i_nSpace]));
244 if (element_active)
245 {
246 isActiveR[offset_u + stride_u*r_l2g[eN_i]] = 1.0;
247 isActiveDOF[u_l2g[eN_i]] = 1.0;
248 }
249 }//i
250 //
251 //save momentum for time history and velocity for subgrid error
252 //save solution for other models
253 //
254
255 q_r[eN_k] = r;
256 q_u[eN_k] = u;
257
258
259 norm = 1.0e-8;
260 for (int I=0;I<nSpace;I++)
261 norm += grad_u[I]*grad_u[I];
262 norm = sqrt(norm);
263 for(int I=0;I<nSpace;I++)
264 q_n[eN_k_nSpace+I] = grad_u[I]/norm;
265 }
266 }
267 void calculateResidual(arguments_dict& args,
268 bool useExact)
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>& x_ref = args.array<double>("x_ref");
275 xt::pyarray<double>& dV_ref = args.array<double>("dV_ref");
276 xt::pyarray<double>& u_trial_ref = args.array<double>("u_trial_ref");
277 xt::pyarray<double>& u_grad_trial_ref = args.array<double>("u_grad_trial_ref");
278 xt::pyarray<double>& u_test_ref = args.array<double>("u_test_ref");
279 xt::pyarray<double>& u_grad_test_ref = args.array<double>("u_grad_test_ref");
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 int nElements_global = args.scalar<int>("nElements_global");
290 double useMetrics = args.scalar<double>("useMetrics");
291 double epsFactHeaviside = args.scalar<double>("epsFactHeaviside");
292 double epsFactDirac = args.scalar<double>("epsFactDirac");
293 double epsFactDiffusion = args.scalar<double>("epsFactDiffusion");
294 xt::pyarray<int>& u_l2g = args.array<int>("u_l2g");
295 xt::pyarray<int>& r_l2g = args.array<int>("r_l2g");
296 xt::pyarray<double>& elementDiameter = args.array<double>("elementDiameter");
297 xt::pyarray<double>& elementBoundaryDiameter = args.array<double>("elementBoundaryDiameter");
298 xt::pyarray<double>& nodeDiametersArray = args.array<double>("nodeDiametersArray");
299 xt::pyarray<double>& u_dof = args.array<double>("u_dof");
300 xt::pyarray<double>& phi_dof = args.array<double>("phi_dof");
301 xt::pyarray<double>& q_phi = args.array<double>("q_phi");
302 xt::pyarray<double>& q_normal_phi = args.array<double>("q_normal_phi");
303 xt::pyarray<double>& ebqe_phi = args.array<double>("ebqe_phi");
304 xt::pyarray<double>& ebqe_normal_phi = args.array<double>("ebqe_normal_phi");
305 xt::pyarray<double>& q_H = args.array<double>("q_H");
306 xt::pyarray<double>& q_u = args.array<double>("q_u");
307 xt::pyarray<double>& q_n = args.array<double>("q_n");
308 xt::pyarray<double>& ebqe_u = args.array<double>("ebqe_u");
309 xt::pyarray<double>& ebqe_n = args.array<double>("ebqe_n");
310 xt::pyarray<double>& q_r = args.array<double>("q_r");
311 xt::pyarray<double>& q_porosity = args.array<double>("q_porosity");
312 int offset_u = args.scalar<int>("offset_u");
313 int stride_u = args.scalar<int>("stride_u");
314 xt::pyarray<double>& globalResidual = args.array<double>("globalResidual");
315 int nExteriorElementBoundaries_global = args.scalar<int>("nExteriorElementBoundaries_global");
316 xt::pyarray<int>& exteriorElementBoundariesArray = args.array<int>("exteriorElementBoundariesArray");
317 xt::pyarray<int>& elementBoundariesArray = args.array<int>("elementBoundariesArray");
318 xt::pyarray<int>& elementBoundaryElementsArray = args.array<int>("elementBoundaryElementsArray");
319 xt::pyarray<int>& elementBoundaryLocalElementBoundariesArray = args.array<int>("elementBoundaryLocalElementBoundariesArray");
320 xt::pyarray<double>& ebqe_phi_s = args.array<double>("ebqe_phi_s");
321 double ghost_penalty_constant = args.scalar<double>("ghost_penalty_constant");
322 const xt::pyarray<double>& phi_solid = args.array<double>("phi_solid");
323 xt::pyarray<double>& phi_solid_nodes = args.array<double>("phi_solid_nodes");
324 bool useExact_s = args.scalar<int>("useExact_s");
325 xt::pyarray<double>& isActiveR = args.array<double>("isActiveR");
326 xt::pyarray<double>& isActiveDOF = args.array<double>("isActiveDOF");
327 xt::pyarray<int>& isActiveElement = args.array<int>("isActiveElement");
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 gf.useExact = useExact;
339 gf_s.useExact = useExact_s;
340 cutfem_boundaries.clear();
341 cutfem_local_boundaries.clear();
342 for(int eN=0;eN<nElements_global;eN++)
343 {
344 //declare local storage for element residual and initialize
345 double elementResidual_u[nDOF_test_element], element_u[nDOF_trial_element], element_phi[nDOF_trial_element], element_phi_s[nDOF_mesh_trial_element];
346 bool element_active=false;
347 isActiveElement[eN]=0;
348 for (int i=0;i<nDOF_test_element;i++)
349 {
350 int eN_i=eN*nDOF_test_element+i;
351 element_u[i] = u_dof.data()[u_l2g.data()[eN_i]];
352 element_phi[i] = phi_dof.data()[u_l2g.data()[eN_i]] + element_u[i];
353 element_phi_s[i] = phi_solid_nodes.data()[u_l2g.data()[eN_i]];
354 }//i
355 double element_nodes[nDOF_mesh_trial_element*3];
356 for (int i=0;i<nDOF_mesh_trial_element;i++)
357 {
358 int eN_i=eN*nDOF_mesh_trial_element+i;
359 for(int I=0;I<3;I++)
360 element_nodes[i*3 + I] = mesh_dof.data()[mesh_l2g.data()[eN_i]*3 + I];
361 }//i
362 gf.calculate(element_phi, element_nodes, x_ref.data(),false);
363 int icase_s = gf_s.calculate(element_phi_s, element_nodes, x_ref.data(),false);
364 if (icase_s == 0)
365 {
366 element_active=true;
367 isActiveElement[eN]=1;
368 //only works for simplices
369 for (int ebN_element=0;ebN_element < nDOF_mesh_trial_element; ebN_element++)
370 {
371 const int ebN = elementBoundariesArray.data()[eN*nDOF_mesh_trial_element+ebN_element];
372 //internal and actually a cut edge
373 //if (elementBoundaryElementsArray.data()[ebN*2+1] != -1 && (ebN < nElementBoundaries_owned) && element_phi_s[(ebN_element+1)%nDOF_mesh_trial_element]*element_phi_s[(ebN_element+2)%nDOF_mesh_trial_element] < 0.0)
374 if (elementBoundaryElementsArray[ebN*2+1] != -1 && element_phi_s[(ebN_element+1)%nDOF_mesh_trial_element]*element_phi_s[(ebN_element+2)%nDOF_mesh_trial_element] <= 0.0)
375 {
376 cutfem_boundaries.insert(ebN);
377 if (elementBoundaryElementsArray[ebN*2 + 0] == eN)
378 cutfem_local_boundaries[ebN] = ebN_element;
379 }
380 }
381 }
382 else if (icase_s == 1)
383 {
384 element_active=true;
385 isActiveElement[eN]=1;
386 }
387 calculateElementResidual(mesh_trial_ref.data(),
388 mesh_grad_trial_ref.data(),
389 mesh_dof.data(),
390 mesh_l2g.data(),
391 dV_ref.data(),
392 u_trial_ref.data(),
393 u_grad_trial_ref.data(),
394 u_test_ref.data(),
395 u_grad_test_ref.data(),
396 mesh_trial_trace_ref.data(),
397 mesh_grad_trial_trace_ref.data(),
398 dS_ref.data(),
399 u_trial_trace_ref.data(),
400 u_grad_trial_trace_ref.data(),
401 u_test_trace_ref.data(),
402 u_grad_test_trace_ref.data(),
403 normal_ref.data(),
404 boundaryJac_ref.data(),
405 nElements_global,
406 useMetrics,
407 epsFactHeaviside,
408 epsFactDirac,
409 epsFactDiffusion,
410 u_l2g.data(),
411 r_l2g.data(),
412 elementDiameter.data(),
413 nodeDiametersArray.data(),
414 u_dof.data(),
415 q_phi.data(),
416 q_normal_phi.data(),
417 ebqe_phi.data(),
418 ebqe_normal_phi.data(),
419 q_H.data(),
420 q_u.data(),
421 q_n.data(),
422 ebqe_u.data(),
423 ebqe_n.data(),
424 q_r.data(),
425 q_porosity.data(),
426 offset_u,stride_u,
427 elementResidual_u,
428 nExteriorElementBoundaries_global,
429 exteriorElementBoundariesArray.data(),
430 elementBoundaryElementsArray.data(),
431 elementBoundaryLocalElementBoundariesArray.data(),
432 element_u,
433 eN,
434 element_active,
435 isActiveR.data(),
436 isActiveDOF.data(),
437 phi_solid.data());
438 //
439 //load element into global residual and save element residual
440 //
441 for(int i=0;i<nDOF_test_element;i++)
442 {
443 int eN_i=eN*nDOF_test_element+i;
444
445 globalResidual.data()[offset_u+stride_u*r_l2g.data()[eN_i]]+=elementResidual_u[i];
446 }//i
447 }//elements
448 std::set<int>::iterator it=cutfem_boundaries.begin();
449 while(it!=cutfem_boundaries.end())
450 {
451 if(isActiveElement[elementBoundaryElementsArray[(*it)*2+0]] && isActiveElement[elementBoundaryElementsArray[(*it)*2+1]])
452 {
453 std::map<int,double> DW_Dn_jump;
454 double gamma_cutfem=ghost_penalty_constant,
455 h_cutfem=elementBoundaryDiameter.data()[*it];
456 int eN_nDOF_trial_element = elementBoundaryElementsArray.data()[(*it)*2+0]*nDOF_trial_element;
457 //See Massing Schott Wall 2018
458 //double norm_v=0.0;
459 //for (int i_offset=1;i_offset<nDOF_trial_element;i_offset++)//MSW18 is just on face, so trying to just use face dof
460 // {
461 // int i = (cutfem_local_boundaries[*it] + i_offset)%nDOF_trial_element;
462 // double u=u_old_dof.data()[vel_l2g.data()[eN_nDOF_trial_element+i]];
463 // v=v_old_dof.data()[vel_l2g.data()[eN_nDOF_v_trial_element+i]],
464 // w=w_old_dof.data()[vel_l2g.data()[eN_nDOF_v_trial_element+i]];
465 // norm_v=fmax(norm_v,sqrt(u*u+v*v+w*w));
466 // }
467 //double gamma_v_dim = rho_0*(nu_0 + norm_v*h_cutfem + alphaBDF*h_cutfem*h_cutfem);
468 //gamma_cutfem_p *= h_cutfem*h_cutfem/gamma_v_dim;
469 //if (NONCONSERVATIVE_FORM)
470 // gamma_cutfem*=gamma_v_dim;
471 //else
472 // gamma_cutfem*=(gamma_v_dim/rho_0);
473 for (int kb=0;kb<nQuadraturePoints_elementBoundary;kb++)
474 {
475 double Du_Dn_jump=0.0, dS;
476 for (int eN_side=0;eN_side < 2; eN_side++)
477 {
478 int ebN = *it,
479 eN = elementBoundaryElementsArray.data()[ebN*2+eN_side];
480 for (int i=0;i<nDOF_test_element;i++)
481 {
482 DW_Dn_jump[r_l2g.data()[eN*nDOF_test_element+i]] = 0.0;
483 }
484 }
485 for (int eN_side=0;eN_side < 2; eN_side++)
486 {
487 int ebN = *it,
488 eN = elementBoundaryElementsArray[ebN*2+eN_side],
489 ebN_local = elementBoundaryLocalElementBoundariesArray[ebN*2+eN_side],
490 eN_nDOF_trial_element = eN*nDOF_trial_element,
491 ebN_local_kb = ebN_local*nQuadraturePoints_elementBoundary+kb,
492 ebN_local_kb_nSpace = ebN_local_kb*nSpace;
493 double u_int=0.0,
494 grad_u_int[nSpace],
495 jac_int[nSpace*nSpace],
496 jacDet_int,
497 jacInv_int[nSpace*nSpace],
498 boundaryJac[nSpace*(nSpace-1)],
499 metricTensor[(nSpace-1)*(nSpace-1)],
500 metricTensorDetSqrt,
501 u_test_dS[nDOF_test_element],
502 u_grad_trial_trace[nDOF_trial_element*nSpace],
503 u_grad_test_dS[nDOF_trial_element*nSpace],
504 normal[nSpace],x_int,y_int,z_int,xt_int,yt_int,zt_int,integralScaling,
505 G[nSpace*nSpace],G_dd_G,tr_G,h_phi,h_penalty,penalty;
506 for (int I=0; I<nSpace;I++)
507 grad_u_int[I] = 0.0;
508 //compute information about mapping from reference element to physical element
509 ck.calculateMapping_elementBoundary(eN,
510 ebN_local,
511 kb,
512 ebN_local_kb,
513 mesh_dof.data(),
514 mesh_l2g.data(),
515 mesh_trial_trace_ref.data(),
516 mesh_grad_trial_trace_ref.data(),
517 boundaryJac_ref.data(),
518 jac_int,
519 jacDet_int,
520 jacInv_int,
521 boundaryJac,
522 metricTensor,
523 metricTensorDetSqrt,
524 normal_ref.data(),
525 normal,
526 x_int,y_int,z_int);
527 dS = metricTensorDetSqrt*dS_ref.data()[kb];
528 //compute shape and solution information
529 //shape
530 ck.gradTrialFromRef(&u_grad_trial_trace_ref.data()[ebN_local_kb_nSpace*nDOF_trial_element],jacInv_int,u_grad_trial_trace);
531 //solution and gradients
532 ck.valFromDOF(u_dof.data(),&u_l2g.data()[eN_nDOF_trial_element],&u_trial_trace_ref.data()[ebN_local_kb*nDOF_test_element],u_int);
533 ck.gradFromDOF(u_dof.data(),&u_l2g.data()[eN_nDOF_trial_element],u_grad_trial_trace,grad_u_int);
534 for (int I=0;I<nSpace;I++)
535 {
536 Du_Dn_jump += grad_u_int[I]*normal[I];
537 }
538 for (int i=0;i<nDOF_test_element;i++)
539 {
540 int eN_i = eN*nDOF_test_element + i;
541 for (int I=0;I<nSpace;I++)
542 DW_Dn_jump[r_l2g[eN_i]] += u_grad_trial_trace[i*nSpace+I]*normal[I];
543 }
544 }//eN_side
545 for (std::map<int,double>::iterator W_it=DW_Dn_jump.begin(); W_it!=DW_Dn_jump.end(); ++W_it)
546 {
547 int i_global = W_it->first;
548 double DW_Dn_jump_i = W_it->second;
549 globalResidual.data()[offset_u+stride_u*i_global]+=gamma_cutfem*h_cutfem*Du_Dn_jump*DW_Dn_jump_i*dS;
550 }
551 }//kb
552 it++;
553 }
554 else
555 {
556 it = cutfem_boundaries.erase(it);
557 }
558 }//cutfem element boundaries
559 //
560 //loop over exterior element boundaries to calculate levelset gradient
561 //
562 //ebNE is the Exterior element boundary INdex
563 //ebN is the element boundary INdex
564 //eN is the element index
565 for (int ebNE = 0; ebNE < nExteriorElementBoundaries_global; ebNE++)
566 {
567 int ebN = exteriorElementBoundariesArray.data()[ebNE],
568 eN = elementBoundaryElementsArray.data()[ebN*2+0],
569 ebN_local = elementBoundaryLocalElementBoundariesArray.data()[ebN*2+0];
570 //eN_nDOF_trial_element = eN*nDOF_trial_element;
571 //double elementResidual_u[nDOF_test_element];
572 double element_u[nDOF_trial_element];
573 for (int i=0;i<nDOF_test_element;i++)
574 {
575 int eN_i=eN*nDOF_test_element+i;
576 element_u[i] = u_dof.data()[u_l2g.data()[eN_i]];
577 }//i
578 for (int kb=0;kb<nQuadraturePoints_elementBoundary;kb++)
579 {
580 int ebNE_kb = ebNE*nQuadraturePoints_elementBoundary+kb,
581 ebNE_kb_nSpace = ebNE_kb*nSpace,
582 ebN_local_kb = ebN_local*nQuadraturePoints_elementBoundary+kb,
583 ebN_local_kb_nSpace = ebN_local_kb*nSpace;
584 double u_ext=0.0,
585 grad_u_ext[nSpace],
586 //m_ext=0.0,
587 //dm_ext=0.0,
588 //H_ext=0.0,
589 //dH_ext[nSpace],
590 //flux_ext=0.0,
591 //bc_u_ext=0.0,
592 //bc_grad_u_ext[nSpace],
593 //bc_m_ext=0.0,
594 //bc_dm_ext=0.0,
595 //bc_H_ext=0.0,
596 //bc_dH_ext[nSpace],
597 jac_ext[nSpace*nSpace],
598 jacDet_ext,
599 jacInv_ext[nSpace*nSpace],
600 boundaryJac[nSpace*(nSpace-1)],
601 metricTensor[(nSpace-1)*(nSpace-1)],
602 metricTensorDetSqrt,
603 dS,
604 //u_test_dS[nDOF_test_element],
605 u_grad_trial_trace[nDOF_trial_element*nSpace],
606 normal[nSpace],x_ext,y_ext,z_ext,
607 G[nSpace*nSpace],G_dd_G,tr_G,norm;
608 //
609 //calculate the solution and gradients at quadrature points
610 //
611 ck.calculateMapping_elementBoundary(eN,
612 ebN_local,
613 kb,
614 ebN_local_kb,
615 mesh_dof.data(),
616 mesh_l2g.data(),
617 mesh_trial_trace_ref.data(),
618 mesh_grad_trial_trace_ref.data(),
619 boundaryJac_ref.data(),
620 jac_ext,
621 jacDet_ext,
622 jacInv_ext,
623 boundaryJac,
624 metricTensor,
625 metricTensorDetSqrt,
626 normal_ref.data(),
627 normal,
628 x_ext,y_ext,z_ext);
629 dS = metricTensorDetSqrt*dS_ref.data()[kb];
630 //get the metric tensor
631 //cek todo use symmetry
632 ck.calculateG(jacInv_ext,G,G_dd_G,tr_G);
633 //compute shape and solution information
634 //shape
635 ck.gradTrialFromRef(&u_grad_trial_trace_ref.data()[ebN_local_kb_nSpace*nDOF_trial_element],jacInv_ext,u_grad_trial_trace);
636 //solution and gradients
637 ck.valFromElementDOF(element_u,&u_trial_trace_ref.data()[ebN_local_kb*nDOF_test_element],u_ext);
638 ck.gradFromElementDOF(element_u,u_grad_trial_trace,grad_u_ext);
639
640 ebqe_u.data()[ebNE_kb] = u_ext;
641 norm = 1.0e-8;
642 for (int I=0;I<nSpace;I++)
643 norm += grad_u_ext[I]*grad_u_ext[I];
644 norm = sqrt(norm);
645 for (int I=0;I<nSpace;I++)
646 ebqe_n.data()[ebNE_kb_nSpace+I] = grad_u_ext[I]/norm;
647 }//kb
648 }//ebNE
649 }
650
651 inline void calculateElementJacobian(//element
652 double* mesh_trial_ref,
653 double* mesh_grad_trial_ref,
654 double* mesh_dof,
655 int* mesh_l2g,
656 double* dV_ref,
657 double* u_trial_ref,
658 double* u_grad_trial_ref,
659 double* u_test_ref,
660 double* u_grad_test_ref,
661 //element boundary
662 double* mesh_trial_trace_ref,
663 double* mesh_grad_trial_trace_ref,
664 double* dS_ref,
665 double* u_trial_trace_ref,
666 double* u_grad_trial_trace_ref,
667 double* u_test_trace_ref,
668 double* u_grad_test_trace_ref,
669 double* normal_ref,
670 double* boundaryJac_ref,
671 //physics
672 int nElements_global,
673 double useMetrics,
674 double epsFactHeaviside,
675 double epsFactDirac,
676 double epsFactDiffusion,
677 int* u_l2g,
678 double* elementDiameter,
679 double* nodeDiametersArray,
680 double* u_dof,
681 // double* u_trial,
682 // double* u_grad_trial,
683 // double* u_test_dV,
684 // double* u_grad_test_dV,
685 double* q_phi,
686 double* q_normal_phi,
687 double* q_H,
688 double* q_porosity,
689 double* elementJacobian_u_u,
690 double* element_u,
691 const double* phi_solid,
692 int eN)
693 {
694 for (int i=0;i<nDOF_test_element;i++)
695 for (int j=0;j<nDOF_trial_element;j++)
696 {
697 elementJacobian_u_u[i*nDOF_trial_element+j]=0.0;
698 }
699 double epsHeaviside,epsDirac,epsDiffusion;
700 for (int k=0;k<nQuadraturePoints_element;k++)
701 {
702 int eN_k = eN*nQuadraturePoints_element+k, //index to a scalar at a quadrature point
703 eN_k_nSpace = eN_k*nSpace;
704 //eN_nDOF_trial_element = eN*nDOF_trial_element; //index to a vector at a quadrature point
705 gf.set_quad(k);
706 gf_s.set_quad(k);
707 //declare local storage
708 double u=0.0,
709 grad_u[nSpace],
710 r=0.0,dr=0.0,
711 jac[nSpace*nSpace],
712 jacDet,
713 jacInv[nSpace*nSpace],
714 u_grad_trial[nDOF_trial_element*nSpace],
715 dV,
716 u_test_dV[nDOF_test_element],
717 u_grad_test_dV[nDOF_test_element*nSpace],
718 x,y,z,
719 G[nSpace*nSpace],G_dd_G,tr_G,h_phi;
720 //
721 //calculate solution and gradients at quadrature points
722 //
723 ck.calculateMapping_element(eN,
724 k,
725 mesh_dof,
726 mesh_l2g,
727 mesh_trial_ref,
728 mesh_grad_trial_ref,
729 jac,
730 jacDet,
731 jacInv,
732 x,y,z);
733 ck.calculateH_element(eN,
734 k,
735 nodeDiametersArray,
736 mesh_l2g,
737 mesh_trial_ref,
738 h_phi);
739 //get the physical integration weight
740 dV = fabs(jacDet)*dV_ref[k];
741 ck.calculateG(jacInv,G,G_dd_G,tr_G);
742
743 /* double dir[nSpace]; */
744 /* double norm = 1.0e-8; */
745 /* for (int I=0;I<nSpace;I++) */
746 /* norm += q_normal_phi[eN_k_nSpace+I]*q_normal_phi[eN_k_nSpace+I]; */
747 /* norm = sqrt(norm); */
748 /* for (int I=0;I<nSpace;I++) */
749 /* dir[I] = q_normal_phi[eN_k_nSpace+I]/norm; */
750 /* ck.calculateGScale(G,dir,h_phi); */
751
752
753 //get the trial function gradients
754 ck.gradTrialFromRef(&u_grad_trial_ref[k*nDOF_trial_element*nSpace],jacInv,u_grad_trial);
755 //get the solution
756 ck.valFromElementDOF(element_u,&u_trial_ref[k*nDOF_trial_element],u);
757 //get the solution gradients
758 ck.gradFromElementDOF(element_u,u_grad_trial,grad_u);
759 //precalculate test function products with integration weights
760 for (int j=0;j<nDOF_trial_element;j++)
761 {
762 u_test_dV[j] = u_test_ref[k*nDOF_trial_element+j]*dV;
763 for (int I=0;I<nSpace;I++)
764 {
765 u_grad_test_dV[j*nSpace+I] = u_grad_trial[j*nSpace+I]*dV;//cek warning won't work for Petrov-Galerkin
766 }
767 }
768 //
769 //calculate pde coefficients and derivatives at quadrature points
770 //
771 epsHeaviside=epsFactHeaviside*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter[eN]);
772 epsDirac =epsFactDirac* (useMetrics*h_phi+(1.0-useMetrics)*elementDiameter[eN]);
773 epsDiffusion=epsFactDiffusion*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter[eN]);
774 // *(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter[eN]);
775 const double H_s = gf_s.H(0.0, phi_solid[eN_k]);
776 evaluateCoefficients(epsHeaviside,
777 epsDirac,
778 q_phi[eN_k],
779 q_H[eN_k],
780 u,
781 q_porosity[eN_k],
782 r,
783 dr);
784 for(int i=0;i<nDOF_test_element;i++)
785 {
786 //int eN_k_i=eN_k*nDOF_test_element+i;
787 //int eN_k_i_nSpace=eN_k_i*nSpace;
788 int i_nSpace=i*nSpace;
789 for(int j=0;j<nDOF_trial_element;j++)
790 {
791 //int eN_k_j=eN_k*nDOF_trial_element+j;
792 //int eN_k_j_nSpace = eN_k_j*nSpace;
793 int j_nSpace = j*nSpace;
794 elementJacobian_u_u[i*nDOF_trial_element+j] +=
795 H_s*(ck.ReactionJacobian_weak(dr,u_trial_ref[k*nDOF_trial_element+j],u_test_dV[i]) +
796 ck.NumericalDiffusionJacobian(epsDiffusion,&u_grad_trial[j_nSpace],&u_grad_test_dV[i_nSpace]));
797 }//j
798 }//i
799 }//k
800 }
801
802 void calculateJacobian(arguments_dict& args,
803 bool useExact)
804 {
805 xt::pyarray<double>& mesh_trial_ref = args.array<double>("mesh_trial_ref");
806 xt::pyarray<double>& mesh_grad_trial_ref = args.array<double>("mesh_grad_trial_ref");
807 xt::pyarray<double>& mesh_dof = args.array<double>("mesh_dof");
808 xt::pyarray<int>& mesh_l2g = args.array<int>("mesh_l2g");
809 xt::pyarray<double>& x_ref = args.array<double>("x_ref");
810 xt::pyarray<double>& dV_ref = args.array<double>("dV_ref");
811 xt::pyarray<double>& u_trial_ref = args.array<double>("u_trial_ref");
812 xt::pyarray<double>& u_grad_trial_ref = args.array<double>("u_grad_trial_ref");
813 xt::pyarray<double>& u_test_ref = args.array<double>("u_test_ref");
814 xt::pyarray<double>& u_grad_test_ref = args.array<double>("u_grad_test_ref");
815 xt::pyarray<double>& mesh_trial_trace_ref = args.array<double>("mesh_trial_trace_ref");
816 xt::pyarray<double>& mesh_grad_trial_trace_ref = args.array<double>("mesh_grad_trial_trace_ref");
817 xt::pyarray<double>& dS_ref = args.array<double>("dS_ref");
818 xt::pyarray<double>& u_trial_trace_ref = args.array<double>("u_trial_trace_ref");
819 xt::pyarray<double>& u_grad_trial_trace_ref = args.array<double>("u_grad_trial_trace_ref");
820 xt::pyarray<double>& u_test_trace_ref = args.array<double>("u_test_trace_ref");
821 xt::pyarray<double>& u_grad_test_trace_ref = args.array<double>("u_grad_test_trace_ref");
822 xt::pyarray<double>& normal_ref = args.array<double>("normal_ref");
823 xt::pyarray<double>& boundaryJac_ref = args.array<double>("boundaryJac_ref");
824 int nElements_global = args.scalar<int>("nElements_global");
825 double useMetrics = args.scalar<double>("useMetrics");
826 double epsFactHeaviside = args.scalar<double>("epsFactHeaviside");
827 double epsFactDirac = args.scalar<double>("epsFactDirac");
828 double epsFactDiffusion = args.scalar<double>("epsFactDiffusion");
829 xt::pyarray<int>& u_l2g = args.array<int>("u_l2g");
830 xt::pyarray<int>& r_l2g = args.array<int>("r_l2g");
831 xt::pyarray<int>& elementBoundariesArray = args.array<int>("elementBoundariesArray");
832 xt::pyarray<int>& elementBoundaryElementsArray = args.array<int>("elementBoundaryElementsArray");
833 xt::pyarray<int>& elementBoundaryLocalElementBoundariesArray = args.array<int>("elementBoundaryLocalElementBoundariesArray");
834 xt::pyarray<double>& elementDiameter = args.array<double>("elementDiameter");
835 xt::pyarray<double>& elementBoundaryDiameter = args.array<double>("elementBoundaryDiameter");
836 xt::pyarray<double>& nodeDiametersArray = args.array<double>("nodeDiametersArray");
837 xt::pyarray<double>& u_dof = args.array<double>("u_dof");
838 xt::pyarray<double>& phi_dof = args.array<double>("phi_dof");
839 xt::pyarray<double>& q_phi = args.array<double>("q_phi");
840 xt::pyarray<double>& q_normal_phi = args.array<double>("q_normal_phi");
841 xt::pyarray<double>& q_H = args.array<double>("q_H");
842 xt::pyarray<double>& q_porosity = args.array<double>("q_porosity");
843 xt::pyarray<int>& csrRowIndeces_u_u = args.array<int>("csrRowIndeces_u_u");
844 xt::pyarray<int>& csrColumnOffsets_u_u = args.array<int>("csrColumnOffsets_u_u");
845 xt::pyarray<int>& csrColumnOffsets_eb_u_u = args.array<int>("csrColumnOffsets_eb_u_u");
846 xt::pyarray<double>& globalJacobian = args.array<double>("globalJacobian");
847 xt::pyarray<double>& ebqe_phi_s = args.array<double>("ebqe_phi_s");
848 const xt::pyarray<double>& phi_solid = args.array<double>("phi_solid");
849 double ghost_penalty_constant = args.scalar<double>("ghost_penalty_constant");
850 xt::pyarray<double>& phi_solid_nodes = args.array<double>("phi_solid_nodes");
851 bool useExact_s = args.scalar<int>("useExact_s");
852 xt::pyarray<double>& isActiveR = args.array<double>("isActiveR");
853 xt::pyarray<double>& isActiveDOF = args.array<double>("isActiveDOF");
854 xt::pyarray<int>& isActiveElement = args.array<int>("isActiveElement");
855 //
856 //loop over elements to compute volume integrals and load them into the element Jacobians and global Jacobian
857 //
858 gf.useExact = useExact;
859 gf_s.useExact = useExact_s;
860 for(int eN=0;eN<nElements_global;eN++)
861 {
862 double elementJacobian_u_u[nDOF_test_element*nDOF_trial_element],element_u[nDOF_trial_element],element_phi[nDOF_trial_element],element_phi_s[nDOF_mesh_trial_element];
863 for (int j=0;j<nDOF_trial_element;j++)
864 {
865 int eN_j = eN*nDOF_trial_element+j;
866 element_u[j] = u_dof.data()[u_l2g.data()[eN_j]];
867 element_phi[j] = phi_dof.data()[u_l2g.data()[eN_j]] + element_u[j];
868 element_phi_s[j] = phi_solid_nodes.data()[u_l2g.data()[eN_j]];
869 }
870 double element_nodes[nDOF_mesh_trial_element*3];
871 for (int i=0;i<nDOF_mesh_trial_element;i++)
872 {
873 int eN_i=eN*nDOF_mesh_trial_element+i;
874 for(int I=0;I<3;I++)
875 element_nodes[i*3 + I] = mesh_dof.data()[mesh_l2g.data()[eN_i]*3 + I];
876 }//i
877 gf.calculate(element_phi, element_nodes, x_ref.data(),false);
878 int icase_s = gf_s.calculate(element_phi_s, element_nodes, x_ref.data(), false);
879 calculateElementJacobian(mesh_trial_ref.data(),
880 mesh_grad_trial_ref.data(),
881 mesh_dof.data(),
882 mesh_l2g.data(),
883 dV_ref.data(),
884 u_trial_ref.data(),
885 u_grad_trial_ref.data(),
886 u_test_ref.data(),
887 u_grad_test_ref.data(),
888 mesh_trial_trace_ref.data(),
889 mesh_grad_trial_trace_ref.data(),
890 dS_ref.data(),
891 u_trial_trace_ref.data(),
892 u_grad_trial_trace_ref.data(),
893 u_test_trace_ref.data(),
894 u_grad_test_trace_ref.data(),
895 normal_ref.data(),
896 boundaryJac_ref.data(),
897 nElements_global,
898 useMetrics,
899 epsFactHeaviside,
900 epsFactDirac,
901 epsFactDiffusion,
902 u_l2g.data(),
903 elementDiameter.data(),
904 nodeDiametersArray.data(),
905 u_dof.data(),
906 q_phi.data(),
907 q_normal_phi.data(),
908 q_H.data(),
909 q_porosity.data(),
910 elementJacobian_u_u,
911 element_u,
912 phi_solid.data(),
913 eN);
914 //
915 //load into element Jacobian into global Jacobian
916 //
917 for (int i=0;i<nDOF_test_element;i++)
918 {
919 int eN_i = eN*nDOF_test_element+i;
920 for (int j=0;j<nDOF_trial_element;j++)
921 {
922 int eN_i_j = eN_i*nDOF_trial_element+j;
923
924 globalJacobian.data()[csrRowIndeces_u_u.data()[eN_i] + csrColumnOffsets_u_u.data()[eN_i_j]] += elementJacobian_u_u[i*nDOF_trial_element+j];
925 }//j
926 }//i
927 }//elements
928 std::set<int>::iterator it=cutfem_boundaries.begin();
929 while(it!=cutfem_boundaries.end())
930 {
931 std::map<int,double> DW_Dn_jump;
932 std::map<std::pair<int, int>, int> u_u_nz;
933 double gamma_cutfem=ghost_penalty_constant,
934 h_cutfem=elementBoundaryDiameter.data()[*it];
935 int eN_nDOF_trial_element = elementBoundaryElementsArray.data()[(*it)*2+0]*nDOF_trial_element;
936 //See Massing Schott Wall 2018
937 //double norm_v=0.0;
938 //for (int i_offset=1;i_offset<nDOF_v_trial_element;i_offset++)//MSW18 is just on face
939 // {
940 // int i = (cutfem_local_boundaries[*it] + i_offset)%nDOF_v_trial_element;//cek hack only works for P1
941 // double u=u_old_dof.data()[vel_l2g.data()[eN_nDOF_v_trial_element+i]],
942 // v=v_old_dof.data()[vel_l2g.data()[eN_nDOF_v_trial_element+i]],
943 // w=w_old_dof.data()[vel_l2g.data()[eN_nDOF_v_trial_element+i]];
944 // norm_v=fmax(norm_v,sqrt(u*u+v*v+w*w));
945 // }
947 //gamma_cutfem_p *= h_cutfem*h_cutfem/gamma_v_dim;
948 //if (NONCONSERVATIVE_FORM)
949 // gamma_cutfem*=gamma_v_dim;
950 //else
951 // gamma_cutfem*=(gamma_v_dim/rho_0);
952 for (int kb=0;kb<nQuadraturePoints_elementBoundary;kb++)
953 {
954 double Du_Dn_jump=0.0, dS;
955 for (int eN_side=0;eN_side < 2; eN_side++)
956 {
957 int ebN = *it,
958 eN = elementBoundaryElementsArray.data()[ebN*2+eN_side];
959 for (int i=0;i<nDOF_test_element;i++)
960 {
961 DW_Dn_jump[r_l2g.data()[eN*nDOF_test_element+i]] = 0.0;
962 }
963 }
964 for (int eN_side=0;eN_side < 2; eN_side++)
965 {
966 int ebN = *it,
967 eN = elementBoundaryElementsArray.data()[ebN*2+eN_side],
968 ebN_local = elementBoundaryLocalElementBoundariesArray.data()[ebN*2+eN_side],
969 eN_nDOF_trial_element = eN*nDOF_trial_element,
970 ebN_local_kb = ebN_local*nQuadraturePoints_elementBoundary+kb,
971 ebN_local_kb_nSpace = ebN_local_kb*nSpace;
972 double u_int=0.0,
973 grad_u_int[nSpace],
974 jac_int[nSpace*nSpace],
975 jacDet_int,
976 jacInv_int[nSpace*nSpace],
977 boundaryJac[nSpace*(nSpace-1)],
978 metricTensor[(nSpace-1)*(nSpace-1)],
979 metricTensorDetSqrt,
980 u_test_dS[nDOF_test_element],
981 u_grad_trial_trace[nDOF_trial_element*nSpace],
982 u_grad_test_dS[nDOF_trial_element*nSpace],
983 normal[nSpace],x_int,y_int,z_int,xt_int,yt_int,zt_int,integralScaling,
984 G[nSpace*nSpace],G_dd_G,tr_G,h_phi,h_penalty,penalty;
985 for (int I=0; I<nSpace;I++)
986 grad_u_int[I] = 0.0;
987 //compute information about mapping from reference element to physical element
988
989 ck.calculateMapping_elementBoundary(eN,
990 ebN_local,
991 kb,
992 ebN_local_kb,
993 mesh_dof.data(),
994 mesh_l2g.data(),
995 mesh_trial_trace_ref.data(),
996 mesh_grad_trial_trace_ref.data(),
997 boundaryJac_ref.data(),
998 jac_int,
999 jacDet_int,
1000 jacInv_int,
1001 boundaryJac,
1002 metricTensor,
1003 metricTensorDetSqrt,
1004 normal_ref.data(),
1005 normal,
1006 x_int,y_int,z_int);
1007 dS = metricTensorDetSqrt*dS_ref.data()[kb];
1008 //compute shape and solution information
1009 //shape
1010 ck.gradTrialFromRef(&u_grad_trial_trace_ref.data()[ebN_local_kb_nSpace*nDOF_trial_element],jacInv_int,u_grad_trial_trace);
1011 for (int i=0;i<nDOF_test_element;i++)
1012 {
1013 int eN_i = eN*nDOF_test_element + i;
1014 for (int I=0;I<nSpace;I++)
1015 DW_Dn_jump[r_l2g.data()[eN_i]] += u_grad_trial_trace[i*nSpace+I]*normal[I];
1016 }
1017 }//eN_side
1018 for (int eN_side=0;eN_side < 2; eN_side++)
1019 {
1020 int ebN = *it,
1021 eN = elementBoundaryElementsArray.data()[ebN*2+eN_side];
1022 for (int i=0;i<nDOF_test_element;i++)
1023 {
1024 int eN_i = eN*nDOF_test_element+i;
1025 for (int eN_side2=0;eN_side2 < 2; eN_side2++)
1026 {
1027 int eN2 = elementBoundaryElementsArray.data()[ebN*2+eN_side2];
1028 for (int j=0;j<nDOF_test_element;j++)
1029 {
1030 int eN_i_j = eN_i*nDOF_test_element + j;
1031 int eN2_j = eN2*nDOF_test_element + j;
1032 int ebN_i_j = ebN*4*nDOF_test_X_trial_element +
1033 eN_side*2*nDOF_test_X_trial_element +
1034 eN_side2*nDOF_test_X_trial_element +
1035 i*nDOF_trial_element +
1036 j;
1037 std::pair<int,int> ij = std::make_pair(u_l2g.data()[eN_i], u_l2g.data()[eN2_j]);
1038 if (u_u_nz.count(ij))
1039 {
1040 assert(u_u_nz[ij] == csrRowIndeces_u_u.data()[eN_i] + csrColumnOffsets_eb_u_u.data()[ebN_i_j]);
1041 }
1042 else
1043 u_u_nz[ij] = csrRowIndeces_u_u.data()[eN_i] + csrColumnOffsets_eb_u_u.data()[ebN_i_j];
1044 }
1045 }
1046 }
1047 }
1048 for (std::map<int,double>::iterator Wi_it=DW_Dn_jump.begin(); Wi_it!=DW_Dn_jump.end(); ++Wi_it)
1049 for (std::map<int,double>::iterator Wj_it=DW_Dn_jump.begin(); Wj_it!=DW_Dn_jump.end(); ++Wj_it)
1050 {
1051 int i_global = Wi_it->first,
1052 j_global = Wj_it->first;
1053 double DW_Dn_jump_i = Wi_it->second,
1054 DW_Dn_jump_j = Wj_it->second;
1055 std::pair<int,int> ij = std::make_pair(i_global, j_global);
1056 globalJacobian.data()[u_u_nz.at(ij)] += gamma_cutfem*h_cutfem*DW_Dn_jump_j*DW_Dn_jump_i*dS;
1057 }//i,j
1058 }//kb
1059 it++;
1060 }//cutfem element boundaries
1061 }//computeJacobian
1062 void elementSolve(arguments_dict& args)
1063 {
1064 xt::pyarray<double>& mesh_trial_ref = args.array<double>("mesh_trial_ref");
1065 xt::pyarray<double>& mesh_grad_trial_ref = args.array<double>("mesh_grad_trial_ref");
1066 xt::pyarray<double>& mesh_dof = args.array<double>("mesh_dof");
1067 xt::pyarray<int>& mesh_l2g = args.array<int>("mesh_l2g");
1068 xt::pyarray<double>& dV_ref = args.array<double>("dV_ref");
1069 xt::pyarray<double>& u_trial_ref = args.array<double>("u_trial_ref");
1070 xt::pyarray<double>& u_grad_trial_ref = args.array<double>("u_grad_trial_ref");
1071 xt::pyarray<double>& u_test_ref = args.array<double>("u_test_ref");
1072 xt::pyarray<double>& u_grad_test_ref = args.array<double>("u_grad_test_ref");
1073 xt::pyarray<double>& mesh_trial_trace_ref = args.array<double>("mesh_trial_trace_ref");
1074 xt::pyarray<double>& mesh_grad_trial_trace_ref = args.array<double>("mesh_grad_trial_trace_ref");
1075 xt::pyarray<double>& dS_ref = args.array<double>("dS_ref");
1076 xt::pyarray<double>& u_trial_trace_ref = args.array<double>("u_trial_trace_ref");
1077 xt::pyarray<double>& u_grad_trial_trace_ref = args.array<double>("u_grad_trial_trace_ref");
1078 xt::pyarray<double>& u_test_trace_ref = args.array<double>("u_test_trace_ref");
1079 xt::pyarray<double>& u_grad_test_trace_ref = args.array<double>("u_grad_test_trace_ref");
1080 xt::pyarray<double>& normal_ref = args.array<double>("normal_ref");
1081 xt::pyarray<double>& boundaryJac_ref = args.array<double>("boundaryJac_ref");
1082 int nElements_global = args.scalar<int>("nElements_global");
1083 double useMetrics = args.scalar<double>("useMetrics");
1084 double epsFactHeaviside = args.scalar<double>("epsFactHeaviside");
1085 double epsFactDirac = args.scalar<double>("epsFactDirac");
1086 double epsFactDiffusion = args.scalar<double>("epsFactDiffusion");
1087 xt::pyarray<int>& u_l2g = args.array<int>("u_l2g");
1088 xt::pyarray<int>& r_l2g = args.array<int>("r_l2g");
1089 xt::pyarray<double>& elementDiameter = args.array<double>("elementDiameter");
1090 xt::pyarray<double>& nodeDiametersArray = args.array<double>("nodeDiametersArray");
1091 xt::pyarray<double>& u_dof = args.array<double>("u_dof");
1092 xt::pyarray<double>& q_phi = args.array<double>("q_phi");
1093 xt::pyarray<double>& q_normal_phi = args.array<double>("q_normal_phi");
1094 xt::pyarray<double>& ebqe_phi = args.array<double>("ebqe_phi");
1095 xt::pyarray<double>& ebqe_normal_phi = args.array<double>("ebqe_normal_phi");
1096 xt::pyarray<double>& q_H = args.array<double>("q_H");
1097 xt::pyarray<double>& q_u = args.array<double>("q_u");
1098 xt::pyarray<double>& q_n = args.array<double>("q_n");
1099 xt::pyarray<double>& ebqe_u = args.array<double>("ebqe_u");
1100 xt::pyarray<double>& ebqe_n = args.array<double>("ebqe_n");
1101 xt::pyarray<double>& q_r = args.array<double>("q_r");
1102 xt::pyarray<double>& q_porosity = args.array<double>("q_porosity");
1103 int offset_u = args.scalar<int>("offset_u");
1104 int stride_u = args.scalar<int>("stride_u");
1105 xt::pyarray<double>& globalResidual = args.array<double>("globalResidual");
1106 int nExteriorElementBoundaries_global = args.scalar<int>("nExteriorElementBoundaries_global");
1107 xt::pyarray<int>& exteriorElementBoundariesArray = args.array<int>("exteriorElementBoundariesArray");
1108 xt::pyarray<int>& elementBoundariesArray = args.array<int>("elementBoundariesArray");
1109 xt::pyarray<int>& elementBoundaryElementsArray = args.array<int>("elementBoundaryElementsArray");
1110 xt::pyarray<int>& elementBoundaryLocalElementBoundariesArray = args.array<int>("elementBoundaryLocalElementBoundariesArray");
1111 xt::pyarray<double>& ebqe_phi_s = args.array<double>("ebqe_phi_s");
1112 const xt::pyarray<double>& phi_solid = args.array<double>("phi_solid");
1113 xt::pyarray<double>& phi_solid_nodes = args.array<double>("phi_solid_nodes");
1114 bool useExact_s = args.scalar<int>("useExact_s");
1115 xt::pyarray<double>& isActiveR = args.array<double>("isActiveR");
1116 xt::pyarray<double>& isActiveDOF = args.array<double>("isActiveDOF");
1117 xt::pyarray<int>& isActiveElement = args.array<int>("isActiveElement");
1118 int maxIts = args.scalar<int>("maxIts");
1119 double atol = args.scalar<double>("atol");
1120 //
1121 //loop over elements to compute volume integrals and load them into element and global residual
1122 //
1123 //eN is the element index
1124 //eN_k is the quadrature point index for a scalar
1125 //eN_k_nSpace is the quadrature point index for a vector
1126 //eN_i is the element test function index
1127 //eN_j is the element trial function index
1128 //eN_k_j is the quadrature point index for a trial function
1129 //eN_k_i is the quadrature point index for a trial function
1130 for(int eN=0;eN<nElements_global;eN++)
1131 {
1132 //declare local storage for element residual and initialize
1133 double element_u[nDOF_test_element],
1134 element_du[nDOF_test_element],
1135 elementResidual_u[nDOF_test_element],
1136 elementJacobian_u_u[nDOF_test_element*nDOF_trial_element],scale=1.0;
1137 PROTEUS_LAPACK_INTEGER elementPivots[nDOF_test_element],
1138 elementColPivots[nDOF_test_element];
1139 //double epsHeaviside,epsDirac,epsDiffusion;
1140 bool element_active=true;
1141 for (int i=0;i<nDOF_test_element;i++)
1142 {
1143 element_u[i]=0.0;
1144 }//i
1145 calculateElementResidual(mesh_trial_ref.data(),
1146 mesh_grad_trial_ref.data(),
1147 mesh_dof.data(),
1148 mesh_l2g.data(),
1149 dV_ref.data(),
1150 u_trial_ref.data(),
1151 u_grad_trial_ref.data(),
1152 u_test_ref.data(),
1153 u_grad_test_ref.data(),
1154 mesh_trial_trace_ref.data(),
1155 mesh_grad_trial_trace_ref.data(),
1156 dS_ref.data(),
1157 u_trial_trace_ref.data(),
1158 u_grad_trial_trace_ref.data(),
1159 u_test_trace_ref.data(),
1160 u_grad_test_trace_ref.data(),
1161 normal_ref.data(),
1162 boundaryJac_ref.data(),
1163 nElements_global,
1164 useMetrics,
1165 epsFactHeaviside,
1166 epsFactDirac,
1167 epsFactDiffusion,
1168 u_l2g.data(),
1169 r_l2g.data(),
1170 elementDiameter.data(),
1171 nodeDiametersArray.data(),
1172 u_dof.data(),
1173 q_phi.data(),
1174 q_normal_phi.data(),
1175 ebqe_phi.data(),
1176 ebqe_normal_phi.data(),
1177 q_H.data(),
1178 q_u.data(),
1179 q_n.data(),
1180 ebqe_u.data(),
1181 ebqe_n.data(),
1182 q_r.data(),
1183 q_porosity.data(),
1184 offset_u,stride_u,
1185 elementResidual_u,
1186 nExteriorElementBoundaries_global,
1187 exteriorElementBoundariesArray.data(),
1188 elementBoundaryElementsArray.data(),
1189 elementBoundaryLocalElementBoundariesArray.data(),
1190 element_u,
1191 eN,
1192 element_active,
1193 isActiveR.data(),
1194 isActiveDOF.data(),
1195 phi_solid.data());
1196 //compute l2 norm
1197 double resNorm=0.0;
1198 for (int i=0;i<nDOF_test_element;i++)
1199 {
1200 resNorm += elementResidual_u[i];
1201 }//i
1202 resNorm = fabs(resNorm);
1203 //now do Newton
1204 int its=0;
1205 //std::cout<<"element "<<eN<<std::endl;
1206 //std::cout<<"resNorm0 "<<resNorm<<std::endl;
1207 while (resNorm >= atol && its < maxIts)
1208 {
1209 its+=1;
1210 calculateElementJacobian(mesh_trial_ref.data(),
1211 mesh_grad_trial_ref.data(),
1212 mesh_dof.data(),
1213 mesh_l2g.data(),
1214 dV_ref.data(),
1215 u_trial_ref.data(),
1216 u_grad_trial_ref.data(),
1217 u_test_ref.data(),
1218 u_grad_test_ref.data(),
1219 mesh_trial_trace_ref.data(),
1220 mesh_grad_trial_trace_ref.data(),
1221 dS_ref.data(),
1222 u_trial_trace_ref.data(),
1223 u_grad_trial_trace_ref.data(),
1224 u_test_trace_ref.data(),
1225 u_grad_test_trace_ref.data(),
1226 normal_ref.data(),
1227 boundaryJac_ref.data(),
1228 nElements_global,
1229 useMetrics,
1230 epsFactHeaviside,
1231 epsFactDirac,
1232 epsFactDiffusion,
1233 u_l2g.data(),
1234 elementDiameter.data(),
1235 nodeDiametersArray.data(),
1236 u_dof.data(),
1237 q_phi.data(),
1238 q_normal_phi.data(),
1239 q_H.data(),
1240 q_porosity.data(),
1241 elementJacobian_u_u,
1242 element_u,
1243 phi_solid.data(),
1244 eN);
1245 for (int i=0;i<nDOF_test_element;i++)
1246 {
1247 element_du[i] = -elementResidual_u[i];
1248 elementPivots[i] = ((PROTEUS_LAPACK_INTEGER)0);
1249 elementColPivots[i]=((PROTEUS_LAPACK_INTEGER)0);
1250 /* std::cout<<"element jacobian"<<std::endl; */
1251 /* for (int j=0;j<nDOF_test_element;j++) */
1252 /* { */
1253 /* std::cout<<elementJacobian_u_u[i*nDOF_trial_element+j]<<'\t'; */
1254 /* } */
1255 /* std::cout<<std::endl; */
1256 }//i
1257 //factor
1258 PROTEUS_LAPACK_INTEGER La_N=((PROTEUS_LAPACK_INTEGER)nDOF_test_element),
1259 INFO=0;
1260 dgetc2_(&La_N,
1261 elementJacobian_u_u,
1262 &La_N,
1263 elementPivots,
1264 elementColPivots,
1265 &INFO);
1266 //solve
1267 dgesc2_(&La_N,
1268 elementJacobian_u_u,
1269 &La_N,
1270 element_du,
1271 elementPivots,
1272 elementColPivots,
1273 &scale);
1274 double resNormNew = resNorm,lambda=1.0;
1275 int lsIts=0;
1276 while (resNormNew > 0.99*resNorm && lsIts < 100)
1277 {
1278 //apply correction
1279 for (int i=0;i<nDOF_test_element;i++)
1280 {
1281 element_u[i] += lambda*element_du[i];
1282 }//i
1283 lambda /= 2.0;
1284 //compute new residual
1285 calculateElementResidual(mesh_trial_ref.data(),
1286 mesh_grad_trial_ref.data(),
1287 mesh_dof.data(),
1288 mesh_l2g.data(),
1289 dV_ref.data(),
1290 u_trial_ref.data(),
1291 u_grad_trial_ref.data(),
1292 u_test_ref.data(),
1293 u_grad_test_ref.data(),
1294 mesh_trial_trace_ref.data(),
1295 mesh_grad_trial_trace_ref.data(),
1296 dS_ref.data(),
1297 u_trial_trace_ref.data(),
1298 u_grad_trial_trace_ref.data(),
1299 u_test_trace_ref.data(),
1300 u_grad_test_trace_ref.data(),
1301 normal_ref.data(),
1302 boundaryJac_ref.data(),
1303 nElements_global,
1304 useMetrics,
1305 epsFactHeaviside,
1306 epsFactDirac,
1307 epsFactDiffusion,
1308 u_l2g.data(),
1309 r_l2g.data(),
1310 elementDiameter.data(),
1311 nodeDiametersArray.data(),
1312 u_dof.data(),
1313 q_phi.data(),
1314 q_normal_phi.data(),
1315 ebqe_phi.data(),
1316 ebqe_normal_phi.data(),
1317 q_H.data(),
1318 q_u.data(),
1319 q_n.data(),
1320 ebqe_u.data(),
1321 ebqe_n.data(),
1322 q_r.data(),
1323 q_porosity.data(),
1324 offset_u,stride_u,
1325 elementResidual_u,
1326 nExteriorElementBoundaries_global,
1327 exteriorElementBoundariesArray.data(),
1328 elementBoundaryElementsArray.data(),
1329 elementBoundaryLocalElementBoundariesArray.data(),
1330 element_u,
1331 eN,
1332 element_active,
1333 isActiveR.data(),
1334 isActiveDOF.data(),
1335 phi_solid.data());
1336
1337 lsIts +=1;
1338 //compute l2 norm
1339 resNormNew=0.0;
1340 for (int i=0;i<nDOF_test_element;i++)
1341 {
1342 resNormNew += elementResidual_u[i];
1343 std::cout<<"element_u["<<i<<"] "<<element_u[i]<<std::endl;
1344 std::cout<<"elementResidual_u["<<i<<"] "<<elementResidual_u[i]<<std::endl;
1345 }//i
1346 resNormNew = fabs(resNormNew);
1347 }
1348 resNorm = resNormNew;
1349 std::cout<<"INFO "<<INFO<<std::endl;
1350 std::cout<<"resNorm["<<its<<"] "<<resNorm<<std::endl;
1351 }
1352 }//elements
1353 }
1354 void elementConstantSolve(arguments_dict& args)
1355 {
1356 xt::pyarray<double>& mesh_trial_ref = args.array<double>("mesh_trial_ref");
1357 xt::pyarray<double>& mesh_grad_trial_ref = args.array<double>("mesh_grad_trial_ref");
1358 xt::pyarray<double>& mesh_dof = args.array<double>("mesh_dof");
1359 xt::pyarray<int>& mesh_l2g = args.array<int>("mesh_l2g");
1360 xt::pyarray<double>& dV_ref = args.array<double>("dV_ref");
1361 xt::pyarray<double>& u_trial_ref = args.array<double>("u_trial_ref");
1362 xt::pyarray<double>& u_grad_trial_ref = args.array<double>("u_grad_trial_ref");
1363 xt::pyarray<double>& u_test_ref = args.array<double>("u_test_ref");
1364 xt::pyarray<double>& u_grad_test_ref = args.array<double>("u_grad_test_ref");
1365 xt::pyarray<double>& mesh_trial_trace_ref = args.array<double>("mesh_trial_trace_ref");
1366 xt::pyarray<double>& mesh_grad_trial_trace_ref = args.array<double>("mesh_grad_trial_trace_ref");
1367 xt::pyarray<double>& dS_ref = args.array<double>("dS_ref");
1368 xt::pyarray<double>& u_trial_trace_ref = args.array<double>("u_trial_trace_ref");
1369 xt::pyarray<double>& u_grad_trial_trace_ref = args.array<double>("u_grad_trial_trace_ref");
1370 xt::pyarray<double>& u_test_trace_ref = args.array<double>("u_test_trace_ref");
1371 xt::pyarray<double>& u_grad_test_trace_ref = args.array<double>("u_grad_test_trace_ref");
1372 xt::pyarray<double>& normal_ref = args.array<double>("normal_ref");
1373 xt::pyarray<double>& boundaryJac_ref = args.array<double>("boundaryJac_ref");
1374 int nElements_global = args.scalar<int>("nElements_global");
1375 double useMetrics = args.scalar<double>("useMetrics");
1376 double epsFactHeaviside = args.scalar<double>("epsFactHeaviside");
1377 double epsFactDirac = args.scalar<double>("epsFactDirac");
1378 double epsFactDiffusion = args.scalar<double>("epsFactDiffusion");
1379 xt::pyarray<int>& u_l2g = args.array<int>("u_l2g");
1380 xt::pyarray<int>& r_l2g = args.array<int>("r_l2g");
1381 xt::pyarray<double>& elementDiameter = args.array<double>("elementDiameter");
1382 xt::pyarray<double>& nodeDiametersArray = args.array<double>("nodeDiametersArray");
1383 xt::pyarray<double>& u_dof = args.array<double>("u_dof");
1384 xt::pyarray<double>& q_phi = args.array<double>("q_phi");
1385 xt::pyarray<double>& q_normal_phi = args.array<double>("q_normal_phi");
1386 xt::pyarray<double>& ebqe_phi = args.array<double>("ebqe_phi");
1387 xt::pyarray<double>& ebqe_normal_phi = args.array<double>("ebqe_normal_phi");
1388 xt::pyarray<double>& q_H = args.array<double>("q_H");
1389 xt::pyarray<double>& q_u = args.array<double>("q_u");
1390 xt::pyarray<double>& q_n = args.array<double>("q_n");
1391 xt::pyarray<double>& ebqe_u = args.array<double>("ebqe_u");
1392 xt::pyarray<double>& ebqe_n = args.array<double>("ebqe_n");
1393 xt::pyarray<double>& q_r = args.array<double>("q_r");
1394 xt::pyarray<double>& q_porosity = args.array<double>("q_porosity");
1395 int offset_u = args.scalar<int>("offset_u");
1396 int stride_u = args.scalar<int>("stride_u");
1397 xt::pyarray<double>& globalResidual = args.array<double>("globalResidual");
1398 int nExteriorElementBoundaries_global = args.scalar<int>("nExteriorElementBoundaries_global");
1399 xt::pyarray<int>& exteriorElementBoundariesArray = args.array<int>("exteriorElementBoundariesArray");
1400 xt::pyarray<int>& elementBoundaryElementsArray = args.array<int>("elementBoundaryElementsArray");
1401 xt::pyarray<int>& elementBoundaryLocalElementBoundariesArray = args.array<int>("elementBoundaryLocalElementBoundariesArray");
1402 const xt::pyarray<double>& phi_solid = args.array<double>("phi_solid");
1403 xt::pyarray<double>& isActiveR = args.array<double>("isActiveR");
1404 xt::pyarray<double>& isActiveDOF = args.array<double>("isActiveDOF");
1405 xt::pyarray<int>& isActiveElement = args.array<int>("isActiveElement");
1406 int maxIts = args.scalar<int>("maxIts");
1407 double atol = args.scalar<double>("atol");
1408 for(int eN=0;eN<nElements_global;eN++)
1409 {
1410 //declare local storage for element residual and initialize
1411 double element_u[nDOF_test_element],elementConstant_u,
1412 elementResidual_u[nDOF_test_element],elementConstantResidual,
1413 elementJacobian_u_u[nDOF_test_element*nDOF_trial_element],elementConstantJacobian,resNorm;
1414 elementConstant_u=0.0;
1415 bool element_active=true;
1416 for (int i=0;i<nDOF_test_element;i++)
1417 {
1418 element_u[i]=elementConstant_u;
1419 }//i
1420 calculateElementResidual(mesh_trial_ref.data(),
1421 mesh_grad_trial_ref.data(),
1422 mesh_dof.data(),
1423 mesh_l2g.data(),
1424 dV_ref.data(),
1425 u_trial_ref.data(),
1426 u_grad_trial_ref.data(),
1427 u_test_ref.data(),
1428 u_grad_test_ref.data(),
1429 mesh_trial_trace_ref.data(),
1430 mesh_grad_trial_trace_ref.data(),
1431 dS_ref.data(),
1432 u_trial_trace_ref.data(),
1433 u_grad_trial_trace_ref.data(),
1434 u_test_trace_ref.data(),
1435 u_grad_test_trace_ref.data(),
1436 normal_ref.data(),
1437 boundaryJac_ref.data(),
1438 nElements_global,
1439 useMetrics,
1440 epsFactHeaviside,
1441 epsFactDirac,
1442 epsFactDiffusion,
1443 u_l2g.data(),
1444 r_l2g.data(),
1445 elementDiameter.data(),
1446 nodeDiametersArray.data(),
1447 u_dof.data(),
1448 q_phi.data(),
1449 q_normal_phi.data(),
1450 ebqe_phi.data(),
1451 ebqe_normal_phi.data(),
1452 q_H.data(),
1453 q_u.data(),
1454 q_n.data(),
1455 ebqe_u.data(),
1456 ebqe_n.data(),
1457 q_r.data(),
1458 q_porosity.data(),
1459 offset_u,stride_u,
1460 elementResidual_u,
1461 nExteriorElementBoundaries_global,
1462 exteriorElementBoundariesArray.data(),
1463 elementBoundaryElementsArray.data(),
1464 elementBoundaryLocalElementBoundariesArray.data(),
1465 element_u,
1466 eN,
1467 element_active,
1468 isActiveR.data(),
1469 isActiveDOF.data(),
1470 phi_solid.data());
1471
1472 //compute l2 norm
1473 elementConstantResidual=0.0;
1474 for (int i=0;i<nDOF_test_element;i++)
1475 {
1476 elementConstantResidual += elementResidual_u[i];
1477 }//i
1478 resNorm = fabs(elementConstantResidual);
1479 //now do Newton
1480 int its=0;
1481 //std::cout<<"element "<<eN<<std::endl;
1482 //std::cout<<"resNorm0 "<<resNorm<<std::endl;
1483 while (resNorm >= atol && its < maxIts)
1484 {
1485 its+=1;
1486 calculateElementJacobian(mesh_trial_ref.data(),
1487 mesh_grad_trial_ref.data(),
1488 mesh_dof.data(),
1489 mesh_l2g.data(),
1490 dV_ref.data(),
1491 u_trial_ref.data(),
1492 u_grad_trial_ref.data(),
1493 u_test_ref.data(),
1494 u_grad_test_ref.data(),
1495 mesh_trial_trace_ref.data(),
1496 mesh_grad_trial_trace_ref.data(),
1497 dS_ref.data(),
1498 u_trial_trace_ref.data(),
1499 u_grad_trial_trace_ref.data(),
1500 u_test_trace_ref.data(),
1501 u_grad_test_trace_ref.data(),
1502 normal_ref.data(),
1503 boundaryJac_ref.data(),
1504 nElements_global,
1505 useMetrics,
1506 epsFactHeaviside,
1507 epsFactDirac,
1508 epsFactDiffusion,
1509 u_l2g.data(),
1510 elementDiameter.data(),
1511 nodeDiametersArray.data(),
1512 u_dof.data(),
1513 q_phi.data(),
1514 q_normal_phi.data(),
1515 q_H.data(),
1516 q_porosity.data(),
1517 elementJacobian_u_u,
1518 element_u,
1519 phi_solid.data(),
1520 eN);
1521 elementConstantJacobian=0.0;
1522 for (int i=0;i<nDOF_test_element;i++)
1523 {
1524 for (int j=0;j<nDOF_test_element;j++)
1525 {
1526 elementConstantJacobian += elementJacobian_u_u[i*nDOF_trial_element+j];
1527 }
1528 }//i
1529 std::cout<<"elementConstantJacobian "<<elementConstantJacobian<<std::endl;
1530 //apply correction
1531 elementConstant_u -= elementConstantResidual/(elementConstantJacobian+1.0e-8);
1532 for (int i=0;i<nDOF_test_element;i++)
1533 {
1534 element_u[i] = elementConstant_u;
1535 }//i
1536 //compute new residual
1537 calculateElementResidual(mesh_trial_ref.data(),
1538 mesh_grad_trial_ref.data(),
1539 mesh_dof.data(),
1540 mesh_l2g.data(),
1541 dV_ref.data(),
1542 u_trial_ref.data(),
1543 u_grad_trial_ref.data(),
1544 u_test_ref.data(),
1545 u_grad_test_ref.data(),
1546 mesh_trial_trace_ref.data(),
1547 mesh_grad_trial_trace_ref.data(),
1548 dS_ref.data(),
1549 u_trial_trace_ref.data(),
1550 u_grad_trial_trace_ref.data(),
1551 u_test_trace_ref.data(),
1552 u_grad_test_trace_ref.data(),
1553 normal_ref.data(),
1554 boundaryJac_ref.data(),
1555 nElements_global,
1556 useMetrics,
1557 epsFactHeaviside,
1558 epsFactDirac,
1559 epsFactDiffusion,
1560 u_l2g.data(),
1561 r_l2g.data(),
1562 elementDiameter.data(),
1563 nodeDiametersArray.data(),
1564 u_dof.data(),
1565 q_phi.data(),
1566 q_normal_phi.data(),
1567 ebqe_phi.data(),
1568 ebqe_normal_phi.data(),
1569 q_H.data(),
1570 q_u.data(),
1571 q_n.data(),
1572 ebqe_u.data(),
1573 ebqe_n.data(),
1574 q_r.data(),
1575 q_porosity.data(),
1576 offset_u,stride_u,
1577 elementResidual_u,
1578 nExteriorElementBoundaries_global,
1579 exteriorElementBoundariesArray.data(),
1580 elementBoundaryElementsArray.data(),
1581 elementBoundaryLocalElementBoundariesArray.data(),
1582 element_u,
1583 eN,
1584 element_active,
1585 isActiveR.data(),
1586 isActiveDOF.data(),
1587 phi_solid.data());
1588
1589 //compute l2 norm
1590 elementConstantResidual=0.0;
1591 for (int i=0;i<nDOF_test_element;i++)
1592 {
1593 elementConstantResidual += elementResidual_u[i];
1594 }//i
1595 resNorm = fabs(elementConstantResidual);
1596 std::cout<<"resNorm["<<its<<"] "<<resNorm<<std::endl;
1597 }
1598 }//elements
1599 }
1600
1601 std::tuple<double, double> globalConstantRJ(arguments_dict& args)
1602 {
1603 xt::pyarray<double>& mesh_trial_ref = args.array<double>("mesh_trial_ref");
1604 xt::pyarray<double>& mesh_grad_trial_ref = args.array<double>("mesh_grad_trial_ref");
1605 xt::pyarray<double>& mesh_dof = args.array<double>("mesh_dof");
1606 xt::pyarray<int>& mesh_l2g = args.array<int>("mesh_l2g");
1607 xt::pyarray<double>& dV_ref = args.array<double>("dV_ref");
1608 xt::pyarray<double>& u_trial_ref = args.array<double>("u_trial_ref");
1609 xt::pyarray<double>& u_grad_trial_ref = args.array<double>("u_grad_trial_ref");
1610 xt::pyarray<double>& u_test_ref = args.array<double>("u_test_ref");
1611 xt::pyarray<double>& u_grad_test_ref = args.array<double>("u_grad_test_ref");
1612 xt::pyarray<double>& mesh_trial_trace_ref = args.array<double>("mesh_trial_trace_ref");
1613 xt::pyarray<double>& mesh_grad_trial_trace_ref = args.array<double>("mesh_grad_trial_trace_ref");
1614 xt::pyarray<double>& dS_ref = args.array<double>("dS_ref");
1615 xt::pyarray<double>& u_trial_trace_ref = args.array<double>("u_trial_trace_ref");
1616 xt::pyarray<double>& u_grad_trial_trace_ref = args.array<double>("u_grad_trial_trace_ref");
1617 xt::pyarray<double>& u_test_trace_ref = args.array<double>("u_test_trace_ref");
1618 xt::pyarray<double>& u_grad_test_trace_ref = args.array<double>("u_grad_test_trace_ref");
1619 xt::pyarray<double>& normal_ref = args.array<double>("normal_ref");
1620 xt::pyarray<double>& boundaryJac_ref = args.array<double>("boundaryJac_ref");
1621 int nElements_owned = args.scalar<int>("nElements_owned");
1622 double useMetrics = args.scalar<double>("useMetrics");
1623 double epsFactHeaviside = args.scalar<double>("epsFactHeaviside");
1624 double epsFactDirac = args.scalar<double>("epsFactDirac");
1625 double epsFactDiffusion = args.scalar<double>("epsFactDiffusion");
1626 xt::pyarray<int>& u_l2g = args.array<int>("u_l2g");
1627 xt::pyarray<int>& r_l2g = args.array<int>("r_l2g");
1628 xt::pyarray<double>& elementDiameter = args.array<double>("elementDiameter");
1629 xt::pyarray<double>& nodeDiametersArray = args.array<double>("nodeDiametersArray");
1630 xt::pyarray<double>& u_dof = args.array<double>("u_dof");
1631 xt::pyarray<double>& q_phi = args.array<double>("q_phi");
1632 xt::pyarray<double>& q_normal_phi = args.array<double>("q_normal_phi");
1633 xt::pyarray<double>& ebqe_phi = args.array<double>("ebqe_phi");
1634 xt::pyarray<double>& ebqe_normal_phi = args.array<double>("ebqe_normal_phi");
1635 xt::pyarray<double>& q_H = args.array<double>("q_H");
1636 xt::pyarray<double>& q_u = args.array<double>("q_u");
1637 xt::pyarray<double>& q_n = args.array<double>("q_n");
1638 xt::pyarray<double>& ebqe_u = args.array<double>("ebqe_u");
1639 xt::pyarray<double>& ebqe_n = args.array<double>("ebqe_n");
1640 xt::pyarray<double>& q_r = args.array<double>("q_r");
1641 xt::pyarray<double>& q_porosity = args.array<double>("q_porosity");
1642 int offset_u = args.scalar<int>("offset_u");
1643 int stride_u = args.scalar<int>("stride_u");
1644 xt::pyarray<double>& globalResidual = args.array<double>("globalResidual");
1645 int nExteriorElementBoundaries_global = args.scalar<int>("nExteriorElementBoundaries_global");
1646 xt::pyarray<int>& exteriorElementBoundariesArray = args.array<int>("exteriorElementBoundariesArray");
1647 xt::pyarray<int>& elementBoundaryElementsArray = args.array<int>("elementBoundaryElementsArray");
1648 xt::pyarray<int>& elementBoundaryLocalElementBoundariesArray = args.array<int>("elementBoundaryLocalElementBoundariesArray");
1649 const xt::pyarray<double>& phi_solid = args.array<double>("phi_solid");
1650 xt::pyarray<double>& isActiveR = args.array<double>("isActiveR");
1651 xt::pyarray<double>& isActiveDOF = args.array<double>("isActiveDOF");
1652 xt::pyarray<int>& isActiveElement = args.array<int>("isActiveElement");
1653 int maxIts = args.scalar<int>("maxIts");
1654 double atol = args.scalar<double>("atol");
1655 double constant_u = args.scalar<double>("constant_u");
1656 double element_u[nDOF_test_element],
1657 elementResidual_u[nDOF_test_element],
1658 elementJacobian_u_u[nDOF_test_element*nDOF_trial_element];
1659 double constantResidual = 0.0;
1660 double constantJacobian = 0.0;
1661 for (int i=0;i<nDOF_trial_element;i++)
1662 {
1663 element_u[i]=constant_u;
1664 }//i
1665 //compute residual and Jacobian
1666 for(int eN=0;eN<nElements_owned;eN++)
1667 {
1668 bool element_active=true;
1669 calculateElementResidual(mesh_trial_ref.data(),
1670 mesh_grad_trial_ref.data(),
1671 mesh_dof.data(),
1672 mesh_l2g.data(),
1673 dV_ref.data(),
1674 u_trial_ref.data(),
1675 u_grad_trial_ref.data(),
1676 u_test_ref.data(),
1677 u_grad_test_ref.data(),
1678 mesh_trial_trace_ref.data(),
1679 mesh_grad_trial_trace_ref.data(),
1680 dS_ref.data(),
1681 u_trial_trace_ref.data(),
1682 u_grad_trial_trace_ref.data(),
1683 u_test_trace_ref.data(),
1684 u_grad_test_trace_ref.data(),
1685 normal_ref.data(),
1686 boundaryJac_ref.data(),
1687 nElements_owned,
1688 useMetrics,
1689 epsFactHeaviside,
1690 epsFactDirac,
1691 epsFactDiffusion,
1692 u_l2g.data(),
1693 r_l2g.data(),
1694 elementDiameter.data(),
1695 nodeDiametersArray.data(),
1696 u_dof.data(),
1697 q_phi.data(),
1698 q_normal_phi.data(),
1699 ebqe_phi.data(),
1700 ebqe_normal_phi.data(),
1701 q_H.data(),
1702 q_u.data(),
1703 q_n.data(),
1704 ebqe_u.data(),
1705 ebqe_n.data(),
1706 q_r.data(),
1707 q_porosity.data(),
1708 offset_u,stride_u,
1709 elementResidual_u,
1710 nExteriorElementBoundaries_global,
1711 exteriorElementBoundariesArray.data(),
1712 elementBoundaryElementsArray.data(),
1713 elementBoundaryLocalElementBoundariesArray.data(),
1714 element_u,
1715 eN,
1716 element_active,
1717 isActiveR.data(),
1718 isActiveDOF.data(),
1719 phi_solid.data());
1720
1721 //compute l2 norm
1722 for (int i=0;i<nDOF_test_element;i++)
1723 {
1724 constantResidual += elementResidual_u[i];
1725 }//i
1726 calculateElementJacobian(mesh_trial_ref.data(),
1727 mesh_grad_trial_ref.data(),
1728 mesh_dof.data(),
1729 mesh_l2g.data(),
1730 dV_ref.data(),
1731 u_trial_ref.data(),
1732 u_grad_trial_ref.data(),
1733 u_test_ref.data(),
1734 u_grad_test_ref.data(),
1735 mesh_trial_trace_ref.data(),
1736 mesh_grad_trial_trace_ref.data(),
1737 dS_ref.data(),
1738 u_trial_trace_ref.data(),
1739 u_grad_trial_trace_ref.data(),
1740 u_test_trace_ref.data(),
1741 u_grad_test_trace_ref.data(),
1742 normal_ref.data(),
1743 boundaryJac_ref.data(),
1744 nElements_owned,
1745 useMetrics,
1746 epsFactHeaviside,
1747 epsFactDirac,
1748 epsFactDiffusion,
1749 u_l2g.data(),
1750 elementDiameter.data(),
1751 nodeDiametersArray.data(),
1752 u_dof.data(),
1753 q_phi.data(),
1754 q_normal_phi.data(),
1755 q_H.data(),
1756 q_porosity.data(),
1757 elementJacobian_u_u,
1758 element_u,
1759 phi_solid.data(),
1760 eN);
1761 for (int i=0;i<nDOF_test_element;i++)
1762 {
1763 for (int j=0;j<nDOF_test_element;j++)
1764 {
1765 constantJacobian += elementJacobian_u_u[i*nDOF_trial_element+j];
1766 }
1767 }//i
1768 }
1769 return std::tuple<double, double>(constantResidual, constantJacobian);
1770 }
1771
1772 double calculateMass(arguments_dict& args,
1773 bool useExact)
1774 {
1775 xt::pyarray<double>& mesh_trial_ref = args.array<double>("mesh_trial_ref");
1776 xt::pyarray<double>& mesh_grad_trial_ref = args.array<double>("mesh_grad_trial_ref");
1777 xt::pyarray<double>& mesh_dof = args.array<double>("mesh_dof");
1778 xt::pyarray<int>& mesh_l2g = args.array<int>("mesh_l2g");
1779 xt::pyarray<double>& x_ref = args.array<double>("x_ref");
1780 xt::pyarray<double>& dV_ref = args.array<double>("dV_ref");
1781 xt::pyarray<double>& u_trial_ref = args.array<double>("u_trial_ref");
1782 xt::pyarray<double>& u_grad_trial_ref = args.array<double>("u_grad_trial_ref");
1783 xt::pyarray<double>& u_test_ref = args.array<double>("u_test_ref");
1784 xt::pyarray<double>& u_grad_test_ref = args.array<double>("u_grad_test_ref");
1785 xt::pyarray<double>& mesh_trial_trace_ref = args.array<double>("mesh_trial_trace_ref");
1786 xt::pyarray<double>& mesh_grad_trial_trace_ref = args.array<double>("mesh_grad_trial_trace_ref");
1787 xt::pyarray<double>& dS_ref = args.array<double>("dS_ref");
1788 xt::pyarray<double>& u_trial_trace_ref = args.array<double>("u_trial_trace_ref");
1789 xt::pyarray<double>& u_grad_trial_trace_ref = args.array<double>("u_grad_trial_trace_ref");
1790 xt::pyarray<double>& u_test_trace_ref = args.array<double>("u_test_trace_ref");
1791 xt::pyarray<double>& u_grad_test_trace_ref = args.array<double>("u_grad_test_trace_ref");
1792 xt::pyarray<double>& normal_ref = args.array<double>("normal_ref");
1793 xt::pyarray<double>& boundaryJac_ref = args.array<double>("boundaryJac_ref");
1794 int nElements_owned = args.scalar<int>("nElements_owned");
1795 double useMetrics = args.scalar<double>("useMetrics");
1796 double epsFactHeaviside = args.scalar<double>("epsFactHeaviside");
1797 double epsFactDirac = args.scalar<double>("epsFactDirac");
1798 double epsFactDiffusion = args.scalar<double>("epsFactDiffusion");
1799 xt::pyarray<int>& u_l2g = args.array<int>("u_l2g");
1800 xt::pyarray<double>& elementDiameter = args.array<double>("elementDiameter");
1801 xt::pyarray<double>& nodeDiametersArray = args.array<double>("nodeDiametersArray");
1802 xt::pyarray<double>& u_dof = args.array<double>("u_dof");
1803 xt::pyarray<double>& phi_dof = args.array<double>("phi_dof");
1804 xt::pyarray<double>& q_phi = args.array<double>("q_phi");
1805 xt::pyarray<double>& q_normal_phi = args.array<double>("q_normal_phi");
1806 xt::pyarray<double>& ebqe_phi = args.array<double>("ebqe_phi");
1807 xt::pyarray<double>& ebqe_normal_phi = args.array<double>("ebqe_normal_phi");
1808 xt::pyarray<double>& q_H = args.array<double>("q_H");
1809 xt::pyarray<double>& q_u = args.array<double>("q_u");
1810 xt::pyarray<double>& q_n = args.array<double>("q_n");
1811 xt::pyarray<double>& ebqe_u = args.array<double>("ebqe_u");
1812 xt::pyarray<double>& ebqe_n = args.array<double>("ebqe_n");
1813 xt::pyarray<double>& q_r = args.array<double>("q_r");
1814 xt::pyarray<double>& q_porosity = args.array<double>("q_porosity");
1815 int offset_u = args.scalar<int>("offset_u");
1816 int stride_u = args.scalar<int>("stride_u");
1817 xt::pyarray<double>& globalResidual = args.array<double>("globalResidual");
1818 int nExteriorElementBoundaries_global = args.scalar<int>("nExteriorElementBoundaries_global");
1819 xt::pyarray<int>& exteriorElementBoundariesArray = args.array<int>("exteriorElementBoundariesArray");
1820 xt::pyarray<int>& elementBoundaryElementsArray = args.array<int>("elementBoundaryElementsArray");
1821 xt::pyarray<int>& elementBoundaryLocalElementBoundariesArray = args.array<int>("elementBoundaryLocalElementBoundariesArray");
1822 const xt::pyarray<double>& phi_solid = args.array<double>("phi_solid");
1823 xt::pyarray<double>& phi_solid_nodes = args.array<double>("phi_solid_nodes");
1824 bool useExact_s = args.scalar<int>("useExact_s");
1825 double globalMass = 0.0;
1826 gf.useExact=useExact;
1827 gf_s.useExact = useExact_s;
1828 for(int eN=0;eN<nElements_owned;eN++)
1829 {
1830 double epsHeaviside;
1831 //loop over quadrature points and compute integrands
1832 //declare local storage for element residual and initialize
1833 double element_phi[nDOF_trial_element],element_phi_s[nDOF_trial_element];
1834 for (int i=0;i<nDOF_test_element;i++)
1835 {
1836 int eN_i=eN*nDOF_test_element+i;
1837 element_phi[i] = phi_dof.data()[u_l2g.data()[eN_i]];
1838 element_phi_s[i] = phi_solid_nodes.data()[u_l2g.data()[eN_i]];
1839 }//i
1840 double element_nodes[nDOF_mesh_trial_element*3];
1841 for (int i=0;i<nDOF_mesh_trial_element;i++)
1842 {
1843 int eN_i=eN*nDOF_mesh_trial_element+i;
1844 for(int I=0;I<3;I++)
1845 element_nodes[i*3 + I] = mesh_dof.data()[mesh_l2g.data()[eN_i]*3 + I];
1846 }//i
1847 gf.calculate(element_phi, element_nodes, x_ref.data(),false);
1848 int icase_s = gf_s.calculate(element_phi_s, element_nodes, x_ref.data(),false);
1849 for (int k=0;k<nQuadraturePoints_element;k++)
1850 {
1851 //compute indeces and declare local storage
1852 int eN_k = eN*nQuadraturePoints_element+k,
1853 eN_k_nSpace = eN_k*nSpace;
1854 //eN_nDOF_trial_element = eN*nDOF_trial_element;
1855 //double u=0.0,grad_u[nSpace],r=0.0,dr=0.0;
1856 double jac[nSpace*nSpace],
1857 jacDet,
1858 jacInv[nSpace*nSpace],
1859 //u_grad_trial[nDOF_trial_element*nSpace],
1860 //u_test_dV[nDOF_trial_element],
1861 //u_grad_test_dV[nDOF_test_element*nSpace],
1862 dV,x,y,z,
1863 G[nSpace*nSpace],G_dd_G,tr_G,h_phi;
1864 gf.set_quad(k);
1865 gf_s.set_quad(k);
1866 //
1867 //compute solution and gradients at quadrature points
1868 //
1869 ck.calculateMapping_element(eN,
1870 k,
1871 mesh_dof.data(),
1872 mesh_l2g.data(),
1873 mesh_trial_ref.data(),
1874 mesh_grad_trial_ref.data(),
1875 jac,
1876 jacDet,
1877 jacInv,
1878 x,y,z);
1879 ck.calculateH_element(eN,
1880 k,
1881 nodeDiametersArray.data(),
1882 mesh_l2g.data(),
1883 mesh_trial_ref.data(),
1884 h_phi);
1885 //get the physical integration weight
1886 dV = fabs(jacDet)*dV_ref.data()[k];
1887 ck.calculateG(jacInv,G,G_dd_G,tr_G);
1888 /* double dir[nSpace]; */
1889 /* double norm = 1.0e-8; */
1890 /* for (int I=0;I<nSpace;I++) */
1891 /* norm += q_normal_phi.data()[eN_k_nSpace+I]*q_normal_phi.data()[eN_k_nSpace+I]; */
1892 /* norm = sqrt(norm); */
1893 /* for (int I=0;I<nSpace;I++) */
1894 /* dir[I] = q_normal_phi.data()[eN_k_nSpace+I]/norm; */
1895 /* ck.calculateGScale(G,dir,h_phi); */
1896 epsHeaviside=epsFactHeaviside*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter.data()[eN]);
1897 globalMass += gf_s.H(epsHeaviside,phi_solid.data()[eN_k])*q_porosity[eN_k]*gf.H(epsHeaviside,q_phi.data()[eN_k])*dV;
1898 }//k
1899 }//elements
1900 return globalMass;
1901 }
1902
1903 void setMassQuadrature(arguments_dict& args,
1904 bool useExact)
1905 {
1906 xt::pyarray<double>& mesh_trial_ref = args.array<double>("mesh_trial_ref");
1907 xt::pyarray<double>& mesh_grad_trial_ref = args.array<double>("mesh_grad_trial_ref");
1908 xt::pyarray<double>& mesh_dof = args.array<double>("mesh_dof");
1909 xt::pyarray<int>& mesh_l2g = args.array<int>("mesh_l2g");
1910 xt::pyarray<double>& x_ref = args.array<double>("x_ref");
1911 xt::pyarray<double>& dV_ref = args.array<double>("dV_ref");
1912 xt::pyarray<double>& u_trial_ref = args.array<double>("u_trial_ref");
1913 xt::pyarray<double>& u_grad_trial_ref = args.array<double>("u_grad_trial_ref");
1914 xt::pyarray<double>& u_test_ref = args.array<double>("u_test_ref");
1915 xt::pyarray<double>& u_grad_test_ref = args.array<double>("u_grad_test_ref");
1916 xt::pyarray<double>& mesh_trial_trace_ref = args.array<double>("mesh_trial_trace_ref");
1917 xt::pyarray<double>& mesh_grad_trial_trace_ref = args.array<double>("mesh_grad_trial_trace_ref");
1918 xt::pyarray<double>& dS_ref = args.array<double>("dS_ref");
1919 xt::pyarray<double>& u_trial_trace_ref = args.array<double>("u_trial_trace_ref");
1920 xt::pyarray<double>& u_grad_trial_trace_ref = args.array<double>("u_grad_trial_trace_ref");
1921 xt::pyarray<double>& u_test_trace_ref = args.array<double>("u_test_trace_ref");
1922 xt::pyarray<double>& u_grad_test_trace_ref = args.array<double>("u_grad_test_trace_ref");
1923 xt::pyarray<double>& normal_ref = args.array<double>("normal_ref");
1924 xt::pyarray<double>& boundaryJac_ref = args.array<double>("boundaryJac_ref");
1925 int nElements_global = args.scalar<int>("nElements_global");
1926 double useMetrics = args.scalar<double>("useMetrics");
1927 double epsFactHeaviside = args.scalar<double>("epsFactHeaviside");
1928 double epsFactDirac = args.scalar<double>("epsFactDirac");
1929 double epsFactDiffusion = args.scalar<double>("epsFactDiffusion");
1930 xt::pyarray<int>& phi_l2g = args.array<int>("phi_l2g");
1931 xt::pyarray<double>& elementDiameter = args.array<double>("elementDiameter");
1932 xt::pyarray<double>& nodeDiametersArray = args.array<double>("nodeDiametersArray");
1933 xt::pyarray<double>& phi_dof = args.array<double>("phi_dof");
1934 xt::pyarray<double>& q_phi = args.array<double>("q_phi");
1935 xt::pyarray<double>& q_normal_phi = args.array<double>("q_normal_phi");
1936 xt::pyarray<double>& ebqe_phi = args.array<double>("ebqe_phi");
1937 xt::pyarray<double>& ebqe_normal_phi = args.array<double>("ebqe_normal_phi");
1938 xt::pyarray<double>& q_H = args.array<double>("q_H");
1939 xt::pyarray<double>& q_u = args.array<double>("q_u");
1940 xt::pyarray<double>& q_n = args.array<double>("q_n");
1941 xt::pyarray<double>& ebqe_u = args.array<double>("ebqe_u");
1942 xt::pyarray<double>& ebqe_n = args.array<double>("ebqe_n");
1943 xt::pyarray<double>& q_r = args.array<double>("q_r");
1944 xt::pyarray<double>& q_porosity = args.array<double>("q_porosity");
1945 int offset_u = args.scalar<int>("offset_u");
1946 int stride_u = args.scalar<int>("stride_u");
1947 xt::pyarray<double>& globalResidual = args.array<double>("globalResidual");
1948 int nExteriorElementBoundaries_global = args.scalar<int>("nExteriorElementBoundaries_global");
1949 xt::pyarray<int>& exteriorElementBoundariesArray = args.array<int>("exteriorElementBoundariesArray");
1950 xt::pyarray<int>& elementBoundaryElementsArray = args.array<int>("elementBoundaryElementsArray");
1951 xt::pyarray<int>& elementBoundaryLocalElementBoundariesArray = args.array<int>("elementBoundaryLocalElementBoundariesArray");
1952 xt::pyarray<double>& H_dof = args.array<double>("H_dof");
1953 gf.useExact=useExact;
1954 gf_nodes.useExact=useExact;
1955 for(int eN=0;eN<nElements_global;eN++)
1956 {
1957 double epsHeaviside;
1958 //loop over quadrature points and compute integrands
1959 //declare local storage for element residual and initialize
1960 double element_phi[nDOF_trial_element];
1961 for (int i=0;i<nDOF_test_element;i++)
1962 {
1963 int eN_i=eN*nDOF_test_element+i;
1964 element_phi[i] = phi_dof.data()[phi_l2g.data()[eN_i]];
1965 }//i
1966 double element_nodes[nDOF_mesh_trial_element*3];
1967 for (int i=0;i<nDOF_mesh_trial_element;i++)
1968 {
1969 int eN_i=eN*nDOF_mesh_trial_element+i;
1970 for(int I=0;I<3;I++)
1971 element_nodes[i*3 + I] = mesh_dof.data()[mesh_l2g.data()[eN_i]*3 + I];
1972 }//i
1973 gf.calculate(element_phi, element_nodes, x_ref.data(),false);
1974 gf_nodes.calculate(element_phi, element_nodes, element_nodes,false);
1975 for (int k=0;k<nQuadraturePoints_element;k++)
1976 {
1977 //compute indeces and declare local storage
1978 int eN_k = eN*nQuadraturePoints_element+k,
1979 eN_k_nSpace = eN_k*nSpace;
1980 //eN_nDOF_trial_element = eN*nDOF_trial_element;
1981 //double u=0.0,grad_u[nSpace],r=0.0,dr=0.0;
1982 double jac[nSpace*nSpace],
1983 jacDet,
1984 jacInv[nSpace*nSpace],
1985 //u_grad_trial[nDOF_trial_element*nSpace],
1986 //u_test_dV[nDOF_trial_element],
1987 //u_grad_test_dV[nDOF_test_element*nSpace],
1988 dV,x,y,z,
1989 G[nSpace*nSpace],G_dd_G,tr_G,h_phi;
1990 //
1991 //compute solution and gradients at quadrature points
1992 //
1993 gf.set_quad(k);
1994 ck.calculateMapping_element(eN,
1995 k,
1996 mesh_dof.data(),
1997 mesh_l2g.data(),
1998 mesh_trial_ref.data(),
1999 mesh_grad_trial_ref.data(),
2000 jac,
2001 jacDet,
2002 jacInv,
2003 x,y,z);
2004 ck.calculateH_element(eN,
2005 k,
2006 nodeDiametersArray.data(),
2007 mesh_l2g.data(),
2008 mesh_trial_ref.data(),
2009 h_phi);
2010 //get the physical integration weight
2011 dV = fabs(jacDet)*dV_ref.data()[k];
2012 ck.calculateG(jacInv,G,G_dd_G,tr_G);
2013
2014 /* double dir[nSpace]; */
2015 /* double norm = 1.0e-8; */
2016 /* for (int I=0;I<nSpace;I++) */
2017 /* norm += q_normal_phi.data()[eN_k_nSpace+I]*q_normal_phi.data()[eN_k_nSpace+I]; */
2018 /* norm = sqrt(norm); */
2019 /* for (int I=0;I<nSpace;I++) */
2020 /* dir[I] = q_normal_phi.data()[eN_k_nSpace+I]/norm; */
2021
2022 /* ck.calculateGScale(G,dir,h_phi); */
2023 epsHeaviside=epsFactHeaviside*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter.data()[eN]);
2024 q_H.data()[eN_k] = gf.H(epsHeaviside,q_phi.data()[eN_k]);
2025 }//k
2026 // distribute rhs for mass correction
2027 for (int i=0;i<nDOF_trial_element;i++)
2028 {
2029 gf_nodes.set_quad(i);
2030 int eN_i = eN*nDOF_trial_element + i;
2031 int gi = phi_l2g.data()[eN_i];
2032 epsHeaviside = epsFactHeaviside*nodeDiametersArray.data()[mesh_l2g.data()[eN_i]];//cek hack, only works if isoparametric, but we can fix by including interpolation points
2033 H_dof.data() [gi] = gf_nodes.H(epsHeaviside,phi_dof.data()[gi]);
2034 }
2035 }//elements
2036 }
2037
2038 void FCTStep(arguments_dict& args)
2039 {
2040 int NNZ = args.scalar<int>("NNZ");
2041 int numDOFs = args.scalar<int>("numDOFs");
2042 xt::pyarray<double>& lumped_mass_matrix = args.array<double>("lumped_mass_matrix");
2043 xt::pyarray<double>& solH = args.array<double>("solH");
2044 xt::pyarray<double>& solL = args.array<double>("solL");
2045 xt::pyarray<double>& limited_solution = args.array<double>("limited_solution");
2046 xt::pyarray<int>& csrRowIndeces_DofLoops = args.array<int>("csrRowIndeces_DofLoops");
2047 xt::pyarray<int>& csrColumnOffsets_DofLoops = args.array<int>("csrColumnOffsets_DofLoops");
2048 xt::pyarray<double>& MassMatrix = args.array<double>("matrix");
2049 Rpos.resize(numDOFs,0.0), Rneg.resize(numDOFs,0.0);
2050 FluxCorrectionMatrix.resize(NNZ,0.0);
2052 // LOOP in DOFs //
2054 int ij=0;
2055 for (int i=0; i<numDOFs; i++)
2056 {
2057 //read some vectors
2058 double solHi = solH.data()[i];
2059 double solLi = solL.data()[i];
2060 double mi = lumped_mass_matrix.data()[i];
2061
2062 double mini=0., maxi=1.0;
2063 double Pposi=0, Pnegi=0;
2064 // LOOP OVER THE SPARSITY PATTERN (j-LOOP)//
2065 for (int offset=csrRowIndeces_DofLoops.data()[i]; offset<csrRowIndeces_DofLoops.data()[i+1]; offset++)
2066 {
2067 int j = csrColumnOffsets_DofLoops.data()[offset];
2068 // i-th row of flux correction matrix
2069 FluxCorrectionMatrix[ij] = ((i==j ? 1. : 0.)*mi - MassMatrix.data()[ij]) * (solH.data()[j]-solHi);
2070
2072 // COMPUTE P VECTORS //
2074 Pposi += FluxCorrectionMatrix[ij]*((FluxCorrectionMatrix[ij] > 0) ? 1. : 0.);
2075 Pnegi += FluxCorrectionMatrix[ij]*((FluxCorrectionMatrix[ij] < 0) ? 1. : 0.);
2076
2077 //update ij
2078 ij+=1;
2079 }
2081 // COMPUTE Q VECTORS //
2083 double Qposi = mi*(maxi-solLi);
2084 double Qnegi = mi*(mini-solLi);
2085
2087 // COMPUTE R VECTORS //
2089 Rpos[i] = ((Pposi==0) ? 1. : std::min(1.0,Qposi/Pposi));
2090 Rneg[i] = ((Pnegi==0) ? 1. : std::min(1.0,Qnegi/Pnegi));
2091 } // i DOFs
2092
2094 // COMPUTE LIMITERS //
2096 ij=0;
2097 for (int i=0; i<numDOFs; i++)
2098 {
2099 double ith_Limiter_times_FluxCorrectionMatrix = 0.;
2100 double Rposi = Rpos[i], Rnegi = Rneg[i];
2101 // LOOP OVER THE SPARSITY PATTERN (j-LOOP)//
2102 for (int offset=csrRowIndeces_DofLoops.data()[i]; offset<csrRowIndeces_DofLoops.data()[i+1]; offset++)
2103 {
2104 int j = csrColumnOffsets_DofLoops.data()[offset];
2105 ith_Limiter_times_FluxCorrectionMatrix +=
2106 ((FluxCorrectionMatrix[ij]>0) ? std::min(Rposi,Rneg[j]) : std::min(Rnegi,Rpos[j]))
2107 * FluxCorrectionMatrix[ij];
2108 //ith_Limiter_times_FluxCorrectionMatrix += FluxCorrectionMatrix[ij];
2109 //update ij
2110 ij+=1;
2111 }
2112 limited_solution.data()[i] = fmax(0.0,solL.data()[i] + 1./lumped_mass_matrix.data()[i]*ith_Limiter_times_FluxCorrectionMatrix);
2113 }
2114 }
2115
2116 // mql. copied from calculateElementJacobian. NOTE: there are some not necessary computations!!!
2117 inline void calculateElementMassMatrix(//element
2118 double* mesh_trial_ref,
2119 double* mesh_grad_trial_ref,
2120 double* mesh_dof,
2121 int* mesh_l2g,
2122 double* dV_ref,
2123 double* u_trial_ref,
2124 double* u_grad_trial_ref,
2125 double* u_test_ref,
2126 double* u_grad_test_ref,
2127 //element boundary
2128 double* mesh_trial_trace_ref,
2129 double* mesh_grad_trial_trace_ref,
2130 double* dS_ref,
2131 double* u_trial_trace_ref,
2132 double* u_grad_trial_trace_ref,
2133 double* u_test_trace_ref,
2134 double* u_grad_test_trace_ref,
2135 double* normal_ref,
2136 double* boundaryJac_ref,
2137 //physics
2138 int nElements_global,
2139 double useMetrics,
2140 double epsFactHeaviside,
2141 double epsFactDirac,
2142 double epsFactDiffusion,
2143 int* u_l2g,
2144 double* elementDiameter,
2145 double* nodeDiametersArray,
2146 double* u_dof,
2147 // double* u_trial,
2148 // double* u_grad_trial,
2149 // double* u_test_dV,
2150 // double* u_grad_test_dV,
2151 double* q_phi,
2152 double* q_normal_phi,
2153 double* q_H,
2154 double* q_porosity,
2155 double* elementMassMatrix,
2156 double* elementLumpedMassMatrix,
2157 double* element_u,
2158 int eN)
2159 {
2160 for (int i=0;i<nDOF_test_element;i++)
2161 {
2162 elementLumpedMassMatrix[i] = 0.0;
2163 for (int j=0;j<nDOF_trial_element;j++)
2164 {
2165 elementMassMatrix[i*nDOF_trial_element+j]=0.0;
2166 }
2167 }
2168 double epsHeaviside,epsDirac,epsDiffusion;
2169 for (int k=0;k<nQuadraturePoints_element;k++)
2170 {
2171 int eN_k = eN*nQuadraturePoints_element+k, //index to a scalar at a quadrature point
2172 eN_k_nSpace = eN_k*nSpace;
2173 //eN_nDOF_trial_element = eN*nDOF_trial_element; //index to a vector at a quadrature point
2174
2175 //declare local storage
2176 double u=0.0,
2177 grad_u[nSpace],
2178 r=0.0,dr=0.0,
2179 jac[nSpace*nSpace],
2180 jacDet,
2181 jacInv[nSpace*nSpace],
2182 u_grad_trial[nDOF_trial_element*nSpace],
2183 dV,
2184 u_test_dV[nDOF_test_element],
2185 u_grad_test_dV[nDOF_test_element*nSpace],
2186 x,y,z,
2187 G[nSpace*nSpace],G_dd_G,tr_G,h_phi;
2188 //
2189 //calculate solution and gradients at quadrature points
2190 //
2191 ck.calculateMapping_element(eN,
2192 k,
2193 mesh_dof,
2194 mesh_l2g,
2195 mesh_trial_ref,
2196 mesh_grad_trial_ref,
2197 jac,
2198 jacDet,
2199 jacInv,
2200 x,y,z);
2201 ck.calculateH_element(eN,
2202 k,
2203 nodeDiametersArray,
2204 mesh_l2g,
2205 mesh_trial_ref,
2206 h_phi);
2207 //get the physical integration weight
2208 dV = fabs(jacDet)*dV_ref[k];
2209 ck.calculateG(jacInv,G,G_dd_G,tr_G);
2210
2211 /* double dir[nSpace]; */
2212 /* double norm = 1.0e-8; */
2213 /* for (int I=0;I<nSpace;I++) */
2214 /* norm += q_normal_phi[eN_k_nSpace+I]*q_normal_phi[eN_k_nSpace+I]; */
2215 /* norm = sqrt(norm); */
2216 /* for (int I=0;I<nSpace;I++) */
2217 /* dir[I] = q_normal_phi[eN_k_nSpace+I]/norm; */
2218 /* ck.calculateGScale(G,dir,h_phi); */
2219
2220
2221 //get the trial function gradients
2222 ck.gradTrialFromRef(&u_grad_trial_ref[k*nDOF_trial_element*nSpace],jacInv,u_grad_trial);
2223 //get the solution
2224 ck.valFromElementDOF(element_u,&u_trial_ref[k*nDOF_trial_element],u);
2225 //get the solution gradients
2226 ck.gradFromElementDOF(element_u,u_grad_trial,grad_u);
2227 //precalculate test function products with integration weights
2228 for (int j=0;j<nDOF_trial_element;j++)
2229 {
2230 u_test_dV[j] = u_test_ref[k*nDOF_trial_element+j]*dV;
2231 for (int I=0;I<nSpace;I++)
2232 {
2233 u_grad_test_dV[j*nSpace+I] = u_grad_trial[j*nSpace+I]*dV;//cek warning won't work for Petrov-Galerkin
2234 }
2235 }
2236 //
2237 //calculate pde coefficients and derivatives at quadrature points
2238 //
2239 epsHeaviside=epsFactHeaviside*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter[eN]);
2240 epsDirac =epsFactDirac* (useMetrics*h_phi+(1.0-useMetrics)*elementDiameter[eN]);
2241 epsDiffusion=epsFactDiffusion*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter[eN]);
2242 // *(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter[eN]);
2243 evaluateCoefficients(epsHeaviside,
2244 epsDirac,
2245 q_phi[eN_k],
2246 q_H[eN_k],
2247 u,
2248 q_porosity[eN_k],
2249 r,
2250 dr);
2251 for(int i=0;i<nDOF_test_element;i++)
2252 {
2253 //int eN_k_i=eN_k*nDOF_test_element+i;
2254 //int eN_k_i_nSpace=eN_k_i*nSpace;
2255 elementLumpedMassMatrix[i] += u_test_dV[i];
2256 for(int j=0;j<nDOF_trial_element;j++)
2257 {
2258 elementMassMatrix[i*nDOF_trial_element+j] += u_trial_ref[k*nDOF_trial_element+j]*u_test_dV[i];
2259 }//j
2260 }//i
2261 }//k
2262 }
2263
2264 void calculateMassMatrix(arguments_dict& args)
2265 {
2266 xt::pyarray<double>& mesh_trial_ref = args.array<double>("mesh_trial_ref");
2267 xt::pyarray<double>& mesh_grad_trial_ref = args.array<double>("mesh_grad_trial_ref");
2268 xt::pyarray<double>& mesh_dof = args.array<double>("mesh_dof");
2269 xt::pyarray<int>& mesh_l2g = args.array<int>("mesh_l2g");
2270 xt::pyarray<double>& dV_ref = args.array<double>("dV_ref");
2271 xt::pyarray<double>& u_trial_ref = args.array<double>("u_trial_ref");
2272 xt::pyarray<double>& u_grad_trial_ref = args.array<double>("u_grad_trial_ref");
2273 xt::pyarray<double>& u_test_ref = args.array<double>("u_test_ref");
2274 xt::pyarray<double>& u_grad_test_ref = args.array<double>("u_grad_test_ref");
2275 xt::pyarray<double>& mesh_trial_trace_ref = args.array<double>("mesh_trial_trace_ref");
2276 xt::pyarray<double>& mesh_grad_trial_trace_ref = args.array<double>("mesh_grad_trial_trace_ref");
2277 xt::pyarray<double>& dS_ref = args.array<double>("dS_ref");
2278 xt::pyarray<double>& u_trial_trace_ref = args.array<double>("u_trial_trace_ref");
2279 xt::pyarray<double>& u_grad_trial_trace_ref = args.array<double>("u_grad_trial_trace_ref");
2280 xt::pyarray<double>& u_test_trace_ref = args.array<double>("u_test_trace_ref");
2281 xt::pyarray<double>& u_grad_test_trace_ref = args.array<double>("u_grad_test_trace_ref");
2282 xt::pyarray<double>& normal_ref = args.array<double>("normal_ref");
2283 xt::pyarray<double>& boundaryJac_ref = args.array<double>("boundaryJac_ref");
2284 int nElements_global = args.scalar<int>("nElements_global");
2285 double useMetrics = args.scalar<double>("useMetrics");
2286 double epsFactHeaviside = args.scalar<double>("epsFactHeaviside");
2287 double epsFactDirac = args.scalar<double>("epsFactDirac");
2288 double epsFactDiffusion = args.scalar<double>("epsFactDiffusion");
2289 xt::pyarray<int>& u_l2g = args.array<int>("u_l2g");
2290 xt::pyarray<double>& elementDiameter = args.array<double>("elementDiameter");
2291 xt::pyarray<double>& nodeDiametersArray = args.array<double>("nodeDiametersArray");
2292 xt::pyarray<double>& u_dof = args.array<double>("u_dof");
2293 xt::pyarray<double>& q_phi = args.array<double>("q_phi");
2294 xt::pyarray<double>& q_normal_phi = args.array<double>("q_normal_phi");
2295 xt::pyarray<double>& q_H = args.array<double>("q_H");
2296 xt::pyarray<double>& q_porosity = args.array<double>("q_porosity");
2297 xt::pyarray<int>& csrRowIndeces_u_u = args.array<int>("csrRowIndeces_u_u");
2298 xt::pyarray<int>& csrColumnOffsets_u_u = args.array<int>("csrColumnOffsets_u_u");
2299 xt::pyarray<double>& globalMassMatrix = args.array<double>("globalMassMatrix");
2300 xt::pyarray<double>& globalLumpedMassMatrix = args.array<double>("globalLumpedMassMatrix");
2301 //
2302 //loop over elements to compute volume integrals and load them into the element Jacobians and global Jacobian
2303 //
2304 for(int eN=0;eN<nElements_global;eN++)
2305 {
2306 double elementMassMatrix[nDOF_test_element*nDOF_trial_element],element_u[nDOF_trial_element], elementLumpedMassMatrix[nDOF_trial_element];
2307 for (int j=0;j<nDOF_trial_element;j++)
2308 {
2309 int eN_j = eN*nDOF_trial_element+j;
2310 element_u[j] = u_dof.data()[u_l2g.data()[eN_j]];
2311 }
2312 calculateElementMassMatrix(mesh_trial_ref.data(),
2313 mesh_grad_trial_ref.data(),
2314 mesh_dof.data(),
2315 mesh_l2g.data(),
2316 dV_ref.data(),
2317 u_trial_ref.data(),
2318 u_grad_trial_ref.data(),
2319 u_test_ref.data(),
2320 u_grad_test_ref.data(),
2321 mesh_trial_trace_ref.data(),
2322 mesh_grad_trial_trace_ref.data(),
2323 dS_ref.data(),
2324 u_trial_trace_ref.data(),
2325 u_grad_trial_trace_ref.data(),
2326 u_test_trace_ref.data(),
2327 u_grad_test_trace_ref.data(),
2328 normal_ref.data(),
2329 boundaryJac_ref.data(),
2330 nElements_global,
2331 useMetrics,
2332 epsFactHeaviside,
2333 epsFactDirac,
2334 epsFactDiffusion,
2335 u_l2g.data(),
2336 elementDiameter.data(),
2337 nodeDiametersArray.data(),
2338 u_dof.data(),
2339 q_phi.data(),
2340 q_normal_phi.data(),
2341 q_H.data(),
2342 q_porosity.data(),
2343 elementMassMatrix,
2344 elementLumpedMassMatrix,
2345 element_u,
2346 eN);
2347 //
2348 //load into element Jacobian into global Jacobian
2349 //
2350 for (int i=0;i<nDOF_test_element;i++)
2351 {
2352 int eN_i = eN*nDOF_test_element+i;
2353 int gi = u_l2g.data()[eN_i];
2354 globalLumpedMassMatrix.data()[gi] += elementLumpedMassMatrix[i];
2355 for (int j=0;j<nDOF_trial_element;j++)
2356 {
2357 int eN_i_j = eN_i*nDOF_trial_element+j;
2358 globalMassMatrix.data()[csrRowIndeces_u_u.data()[eN_i] + csrColumnOffsets_u_u.data()[eN_i_j]] +=
2359 elementMassMatrix[i*nDOF_trial_element+j];
2360 }//j
2361 }//i
2362 }//elements
2363 }//calculate mass matrix
2364
2365 void setMassQuadratureEdgeBasedStabilizationMethods(arguments_dict& args,
2366 bool useExact)
2367 {
2368 xt::pyarray<double>& mesh_trial_ref = args.array<double>("mesh_trial_ref");
2369 xt::pyarray<double>& mesh_grad_trial_ref = args.array<double>("mesh_grad_trial_ref");
2370 xt::pyarray<double>& mesh_dof = args.array<double>("mesh_dof");
2371 xt::pyarray<int>& mesh_l2g = args.array<int>("mesh_l2g");
2372 xt::pyarray<double>& x_ref = args.array<double>("x_ref");
2373 xt::pyarray<double>& dV_ref = args.array<double>("dV_ref");
2374 xt::pyarray<double>& u_trial_ref = args.array<double>("u_trial_ref");
2375 xt::pyarray<double>& u_grad_trial_ref = args.array<double>("u_grad_trial_ref");
2376 xt::pyarray<double>& u_test_ref = args.array<double>("u_test_ref");
2377 xt::pyarray<double>& u_grad_test_ref = args.array<double>("u_grad_test_ref");
2378 xt::pyarray<double>& mesh_trial_trace_ref = args.array<double>("mesh_trial_trace_ref");
2379 xt::pyarray<double>& mesh_grad_trial_trace_ref = args.array<double>("mesh_grad_trial_trace_ref");
2380 xt::pyarray<double>& dS_ref = args.array<double>("dS_ref");
2381 xt::pyarray<double>& u_trial_trace_ref = args.array<double>("u_trial_trace_ref");
2382 xt::pyarray<double>& u_grad_trial_trace_ref = args.array<double>("u_grad_trial_trace_ref");
2383 xt::pyarray<double>& u_test_trace_ref = args.array<double>("u_test_trace_ref");
2384 xt::pyarray<double>& u_grad_test_trace_ref = args.array<double>("u_grad_test_trace_ref");
2385 xt::pyarray<double>& normal_ref = args.array<double>("normal_ref");
2386 xt::pyarray<double>& boundaryJac_ref = args.array<double>("boundaryJac_ref");
2387 int nElements_global = args.scalar<int>("nElements_global");
2388 double useMetrics = args.scalar<double>("useMetrics");
2389 double epsFactHeaviside = args.scalar<double>("epsFactHeaviside");
2390 double epsFactDirac = args.scalar<double>("epsFactDirac");
2391 double epsFactDiffusion = args.scalar<double>("epsFactDiffusion");
2392 xt::pyarray<int>& phi_l2g = args.array<int>("phi_l2g");
2393 xt::pyarray<double>& elementDiameter = args.array<double>("elementDiameter");
2394 xt::pyarray<double>& nodeDiametersArray = args.array<double>("nodeDiametersArray");
2395 xt::pyarray<double>& phi_dof = args.array<double>("phi_dof");
2396 xt::pyarray<double>& q_phi = args.array<double>("q_phi");
2397 xt::pyarray<double>& q_normal_phi = args.array<double>("q_normal_phi");
2398 xt::pyarray<double>& ebqe_phi = args.array<double>("ebqe_phi");
2399 xt::pyarray<double>& ebqe_normal_phi = args.array<double>("ebqe_normal_phi");
2400 xt::pyarray<double>& q_H = args.array<double>("q_H");
2401 xt::pyarray<double>& q_u = args.array<double>("q_u");
2402 xt::pyarray<double>& q_n = args.array<double>("q_n");
2403 xt::pyarray<double>& ebqe_u = args.array<double>("ebqe_u");
2404 xt::pyarray<double>& ebqe_n = args.array<double>("ebqe_n");
2405 xt::pyarray<double>& q_r = args.array<double>("q_r");
2406 xt::pyarray<double>& q_porosity = args.array<double>("q_porosity");
2407 int offset_u = args.scalar<int>("offset_u");
2408 int stride_u = args.scalar<int>("stride_u");
2409 xt::pyarray<double>& globalResidual = args.array<double>("globalResidual");
2410 int nExteriorElementBoundaries_global = args.scalar<int>("nExteriorElementBoundaries_global");
2411 xt::pyarray<int>& exteriorElementBoundariesArray = args.array<int>("exteriorElementBoundariesArray");
2412 xt::pyarray<int>& elementBoundaryElementsArray = args.array<int>("elementBoundaryElementsArray");
2413 xt::pyarray<int>& elementBoundaryLocalElementBoundariesArray = args.array<int>("elementBoundaryLocalElementBoundariesArray");
2414 xt::pyarray<double>& rhs_mass_correction = args.array<double>("rhs_mass_correction");
2415 xt::pyarray<double>& lumped_L2p_vof_mass_correction = args.array<double>("lumped_L2p_vof_mass_correction");
2416 xt::pyarray<double>& lumped_mass_matrix = args.array<double>("lumped_mass_matrix");
2417 int numDOFs = args.scalar<int>("numDOFs");
2418 gf.useExact=useExact;
2419 for(int eN=0;eN<nElements_global;eN++)
2420 {
2421 double element_rhs_mass_correction[nDOF_test_element];
2422 for (int i=0;i<nDOF_test_element;i++)
2423 element_rhs_mass_correction[i] = 0.;
2424 double epsHeaviside;
2425 //loop over quadrature points and compute integrands
2426 double element_phi[nDOF_trial_element];
2427 for (int i=0;i<nDOF_test_element;i++)
2428 {
2429 int eN_i=eN*nDOF_test_element+i;
2430 element_phi[i] = phi_dof.data()[phi_l2g.data()[eN_i]];
2431 }//i
2432 double element_nodes[nDOF_mesh_trial_element*3];
2433 for (int i=0;i<nDOF_mesh_trial_element;i++)
2434 {
2435 int eN_i=eN*nDOF_mesh_trial_element+i;
2436 for(int I=0;I<3;I++)
2437 element_nodes[i*3 + I] = mesh_dof.data()[mesh_l2g.data()[eN_i]*3 + I];
2438 }//i
2439 gf.calculate(element_phi, element_nodes, x_ref.data(),false);
2440 for (int k=0;k<nQuadraturePoints_element;k++)
2441 {
2442 //compute indeces and declare local storage
2443 int eN_k = eN*nQuadraturePoints_element+k,
2444 eN_k_nSpace = eN_k*nSpace;
2445 //eN_nDOF_trial_element = eN*nDOF_trial_element;
2446 //double u=0.0,grad_u[nSpace],r=0.0,dr=0.0;
2447 double jac[nSpace*nSpace],
2448 jacDet,
2449 jacInv[nSpace*nSpace],
2450 //u_grad_trial[nDOF_trial_element*nSpace],
2451 //u_test_dV[nDOF_trial_element],
2452 //u_grad_test_dV[nDOF_test_element*nSpace],
2453 dV,x,y,z,
2454 u_test_dV[nDOF_test_element],
2455 G[nSpace*nSpace],G_dd_G,tr_G,h_phi;
2456 gf.set_quad(k);
2457 //
2458 ck.calculateMapping_element(eN,
2459 k,
2460 mesh_dof.data(),
2461 mesh_l2g.data(),
2462 mesh_trial_ref.data(),
2463 mesh_grad_trial_ref.data(),
2464 jac,
2465 jacDet,
2466 jacInv,
2467 x,y,z);
2468 ck.calculateH_element(eN,
2469 k,
2470 nodeDiametersArray.data(),
2471 mesh_l2g.data(),
2472 mesh_trial_ref.data(),
2473 h_phi);
2474 //get the physical integration weight
2475 dV = fabs(jacDet)*dV_ref.data()[k];
2476 ck.calculateG(jacInv,G,G_dd_G,tr_G);
2477
2478 // precalculate test function times integration weight
2479 for (int j=0;j<nDOF_trial_element;j++)
2480 u_test_dV[j] = u_test_ref.data()[k*nDOF_trial_element+j]*dV;
2481
2482 /* double dir[nSpace]; */
2483 /* double norm = 1.0e-8; */
2484 /* for (int I=0;I<nSpace;I++) */
2485 /* norm += q_normal_phi.data()[eN_k_nSpace+I]*q_normal_phi.data()[eN_k_nSpace+I]; */
2486 /* norm = sqrt(norm); */
2487 /* for (int I=0;I<nSpace;I++) */
2488 /* dir[I] = q_normal_phi.data()[eN_k_nSpace+I]/norm; */
2489
2490 /* ck.calculateGScale(G,dir,h_phi); */
2491 epsHeaviside=epsFactHeaviside*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter.data()[eN]);
2492 q_H.data()[eN_k] = gf.H(epsHeaviside,q_phi.data()[eN_k]);
2493
2494 for (int i=0;i<nDOF_trial_element;i++)
2495 element_rhs_mass_correction [i] += q_porosity.data()[eN_k]*q_H.data()[eN_k]*u_test_dV[i];
2496 }//k
2497 // distribute rhs for mass correction
2498 for (int i=0;i<nDOF_trial_element;i++)
2499 {
2500 int eN_i = eN*nDOF_trial_element + i;
2501 int gi = phi_l2g.data()[eN_i];
2502 rhs_mass_correction.data()[gi] += element_rhs_mass_correction[i];
2503 }
2504 }//elements
2505 // COMPUTE LUMPED L2 PROYJECTION
2506 for (int i=0; i<numDOFs; i++)
2507 {
2508 double mi = lumped_mass_matrix.data()[i];
2509 lumped_L2p_vof_mass_correction.data()[i] = 1./mi*rhs_mass_correction.data()[i];
2510 }
2511 }
2512 };//MCorr
2513
2514 inline MCorr_base* newMCorr(int nSpaceIn,
2515 int nQuadraturePoints_elementIn,
2516 int nDOF_mesh_trial_elementIn,
2517 int nDOF_trial_elementIn,
2518 int nDOF_test_elementIn,
2519 int nQuadraturePoints_elementBoundaryIn,
2520 int CompKernelFlag)
2521 {
2522 if (nSpaceIn == 2)
2523 return proteus::chooseAndAllocateDiscretization2D<MCorr_base,MCorr,CompKernel>(nSpaceIn,
2524 nQuadraturePoints_elementIn,
2525 nDOF_mesh_trial_elementIn,
2526 nDOF_trial_elementIn,
2527 nDOF_test_elementIn,
2528 nQuadraturePoints_elementBoundaryIn,
2529 CompKernelFlag);
2530 else
2531 return proteus::chooseAndAllocateDiscretization<MCorr_base,MCorr,CompKernel>(nSpaceIn,
2532 nQuadraturePoints_elementIn,
2533 nDOF_mesh_trial_elementIn,
2534 nDOF_trial_elementIn,
2535 nDOF_test_elementIn,
2536 nQuadraturePoints_elementBoundaryIn,
2537 CompKernelFlag);
2538 }
2539}//proteus
2540#endif
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
equivalent_polynomials::GeneralizedFunctions_mix::useExact
bool useExact
Definition equivalent_polynomials.h:1289
equivalent_polynomials::GeneralizedFunctions_mix::set_quad
void set_quad(unsigned int q)
Definition equivalent_polynomials.h:1325
equivalent_polynomials::GeneralizedFunctions_mix::H
double H(double eps, double phi)
Definition equivalent_polynomials.h:1337
equivalent_polynomials::GeneralizedFunctions_mix::D
double D(double eps, double phi)
Definition equivalent_polynomials.h:1353
equivalent_polynomials::GeneralizedFunctions_mix::calculate
int calculate(const double *phi_dof, const double *phi_nodes, const double *xi_r, double ma, double mb, double jf, bool isBoundary, bool scale)
Definition equivalent_polynomials.h:1294
proteus::MCorr_base
Definition MCorr.h:33
proteus::MCorr_base::Rpos
std::valarray< double > Rpos
Definition MCorr.h:35
proteus::MCorr_base::globalConstantRJ
virtual std::tuple< double, double > globalConstantRJ(arguments_dict &args)=0
proteus::MCorr_base::calculateMassMatrix
virtual void calculateMassMatrix(arguments_dict &args)=0
proteus::MCorr_base::elementConstantSolve
virtual void elementConstantSolve(arguments_dict &args)=0
proteus::MCorr_base::setMassQuadratureEdgeBasedStabilizationMethods
virtual void setMassQuadratureEdgeBasedStabilizationMethods(arguments_dict &args, bool useExact)=0
proteus::MCorr_base::setMassQuadrature
virtual void setMassQuadrature(arguments_dict &args, bool useExact)=0
proteus::MCorr_base::calculateMass
virtual double calculateMass(arguments_dict &args, bool useExact)=0
proteus::MCorr_base::calculateResidual
virtual void calculateResidual(arguments_dict &args, bool useExact)=0
proteus::MCorr_base::calculateJacobian
virtual void calculateJacobian(arguments_dict &args, bool useExact)=0
proteus::MCorr_base::Rneg
std::valarray< double > Rneg
Definition MCorr.h:35
proteus::MCorr_base::~MCorr_base
virtual ~MCorr_base()
Definition MCorr.h:37
proteus::MCorr_base::FCTStep
virtual void FCTStep(arguments_dict &args)=0
proteus::MCorr_base::FluxCorrectionMatrix
std::valarray< double > FluxCorrectionMatrix
Definition MCorr.h:36
proteus::MCorr_base::elementSolve
virtual void elementSolve(arguments_dict &args)=0
proteus::MCorr::gf_s
GeneralizedFunctions< nSpace, 4, nQuadraturePoints_element, nQuadraturePoints_elementBoundary > gf_s
Definition MCorr.h:67
proteus::MCorr::ck
CompKernelType ck
Definition MCorr.h:64
proteus::MCorr::gf
GeneralizedFunctions< nSpace, 2, nQuadraturePoints_element, nQuadraturePoints_elementBoundary > gf
Definition MCorr.h:65
proteus::MCorr::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, int *r_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_porosity, int offset_u, int stride_u, double *elementResidual_u, int nExteriorElementBoundaries_global, int *exteriorElementBoundariesArray, int *elementBoundaryElementsArray, int *elementBoundaryLocalElementBoundariesArray, double *element_u, int eN, bool element_active, double *isActiveR, double *isActiveDOF, const double *phi_solid)
Definition MCorr.h:87
proteus::MCorr::calculateMassMatrix
void calculateMassMatrix(arguments_dict &args)
Definition MCorr.h:2264
proteus::MCorr::calculateJacobian
void calculateJacobian(arguments_dict &args, bool useExact)
Definition MCorr.h:802
proteus::MCorr::elementConstantSolve
void elementConstantSolve(arguments_dict &args)
Definition MCorr.h:1354
proteus::MCorr::cutfem_local_boundaries
std::map< int, int > cutfem_local_boundaries
Definition MCorr.h:62
proteus::MCorr::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 MCorr.h:74
proteus::MCorr::setMassQuadratureEdgeBasedStabilizationMethods
void setMassQuadratureEdgeBasedStabilizationMethods(arguments_dict &args, bool useExact)
Definition MCorr.h:2365
proteus::MCorr::gf_nodes
GeneralizedFunctions< nSpace, 2, nDOF_trial_element, nQuadraturePoints_elementBoundary > gf_nodes
Definition MCorr.h:66
proteus::MCorr::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_porosity, double *elementJacobian_u_u, double *element_u, const double *phi_solid, int eN)
Definition MCorr.h:651
proteus::MCorr::setMassQuadrature
void setMassQuadrature(arguments_dict &args, bool useExact)
Definition MCorr.h:1903
proteus::MCorr::calculateMass
double calculateMass(arguments_dict &args, bool useExact)
Definition MCorr.h:1772
proteus::MCorr::FCTStep
void FCTStep(arguments_dict &args)
Definition MCorr.h:2038
proteus::MCorr::elementSolve
void elementSolve(arguments_dict &args)
Definition MCorr.h:1062
proteus::MCorr::globalConstantRJ
std::tuple< double, double > globalConstantRJ(arguments_dict &args)
Definition MCorr.h:1601
proteus::MCorr::calculateElementMassMatrix
void calculateElementMassMatrix(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_porosity, double *elementMassMatrix, double *elementLumpedMassMatrix, double *element_u, int eN)
Definition MCorr.h:2117
proteus::MCorr::calculateResidual
void calculateResidual(arguments_dict &args, bool useExact)
Definition MCorr.h:267
proteus::MCorr::nDOF_test_X_trial_element
const int nDOF_test_X_trial_element
Definition MCorr.h:63
proteus::MCorr::cutfem_boundaries
std::set< int > cutfem_boundaries
Definition MCorr.h:61
equivalent_polynomials.h
proteus
Definition ADR.h:19
proteus::GeneralizedFunctions
equivalent_polynomials::GeneralizedFunctions_mix< nSpace, nP_ifem, nP, nQ, nEBQ, true > GeneralizedFunctions
Definition ADR.h:21
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::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::newMCorr
MCorr_base * newMCorr(int nSpaceIn, int nQuadraturePoints_elementIn, int nDOF_mesh_trial_elementIn, int nDOF_trial_elementIn, int nDOF_test_elementIn, int nQuadraturePoints_elementBoundaryIn, int CompKernelFlag)
Definition MCorr.h:2514
dgesc2_
int dgesc2_(int *n, double *a, int *lda, double *rhs, int *ipiv, int *jpiv, double *scale)
dgetc2_
int dgetc2_(int *n, double *a, int *lda, int *ipiv, int *jpiv, int *info)
proteus::arguments_dict
Definition ArgumentsDict.h:70
proteus::arguments_dict::scalar
T & scalar(const std::string &key)
proteus::arguments_dict::array
xt::pyarray< T > & array(const std::string &key)
Generated on for proteus by doxygen 1.18.0