7 #include PROTEUS_LAPACK_H
10 #include "xtensor-python/pyarray.hpp"
12 #define FAST_ASSEMBLY 1
30 template<
class CompKernelType,
32 int nQuadraturePoints_element,
33 int nDOF_mesh_trial_element,
34 int nDOF_trial_element,
35 int nDOF_test_element,
36 int nQuadraturePoints_elementBoundary>
53 H = 0.5*(1.0 +
phi/eps + sin(M_PI*
phi/eps)/M_PI);
62 HI=
phi - eps + 0.5*(eps + 0.5*eps*eps/eps - eps*cos(M_PI*eps/eps)/(M_PI*M_PI)) - 0.5*((-eps) + 0.5*(-eps)*(-eps)/eps - eps*cos(M_PI*(-eps)/eps)/(M_PI*M_PI));
70 HI = 0.5*(
phi + 0.5*
phi*
phi/eps - eps*cos(M_PI*
phi/eps)/(M_PI*M_PI)) - 0.5*((-eps) + 0.5*(-eps)*(-eps)/eps - eps*cos(M_PI*(-eps)/eps)/(M_PI*M_PI));
83 d = 0.5*(1.0 + cos(M_PI*
phi/eps))/eps;
88 const double& epsDirac,
92 const double& porosity,
101 double* mesh_trial_ref,
102 double* mesh_grad_trial_ref,
107 double* u_grad_trial_ref,
109 double* u_grad_test_ref,
111 double* mesh_trial_trace_ref,
112 double* mesh_grad_trial_trace_ref,
114 double* u_trial_trace_ref,
115 double* u_grad_trial_trace_ref,
116 double* u_test_trace_ref,
117 double* u_grad_test_trace_ref,
119 double* boundaryJac_ref,
121 int nElements_global,
123 double epsFactHeaviside,
125 double epsFactDiffusion,
127 double* elementDiameter,
128 double* nodeDiametersArray,
131 double* q_normal_phi,
133 double* ebqe_normal_phi,
141 int offset_u,
int stride_u,
142 double* elementResidual_u,
143 double* elementInterface_lumpedMassMatrix,
144 int nExteriorElementBoundaries_global,
145 int* exteriorElementBoundariesArray,
146 int* elementBoundaryElementsArray,
147 int* elementBoundaryLocalElementBoundariesArray,
151 for (
int i=0;i<nDOF_test_element;i++)
153 elementResidual_u[i]=0.0;
154 elementInterface_lumpedMassMatrix[i]=0.0;
156 double epsHeaviside,epsDirac,epsDiffusion,norm;
158 for (
int k=0;k<nQuadraturePoints_element;k++)
161 register int eN_k = eN*nQuadraturePoints_element+k,
162 eN_k_nSpace = eN_k*nSpace;
164 register double u=0.0,grad_u[nSpace],
168 jacInv[nSpace*nSpace],
169 u_grad_trial[nDOF_trial_element*nSpace],
170 u_test_dV[nDOF_trial_element],
171 u_grad_test_dV[nDOF_test_element*nSpace],
173 G[nSpace*nSpace],G_dd_G,tr_G,h_phi;
177 ck.calculateMapping_element(eN,
187 ck.calculateH_element(eN,
194 dV = fabs(jacDet)*dV_ref[k];
195 ck.calculateG(jacInv,G,G_dd_G,tr_G);
207 ck.gradTrialFromRef(&u_grad_trial_ref[k*nDOF_trial_element*nSpace],jacInv,u_grad_trial);
209 ck.valFromElementDOF(element_u,&u_trial_ref[k*nDOF_trial_element],
u);
211 ck.gradFromElementDOF(element_u,u_grad_trial,grad_u);
213 for (
int j=0;j<nDOF_trial_element;j++)
215 u_test_dV[j] = u_test_ref[k*nDOF_trial_element+j]*dV;
216 for (
int I=0;I<nSpace;I++)
218 u_grad_test_dV[j*nSpace+I] = u_grad_trial[j*nSpace+I]*dV;
224 epsHeaviside = epsFactHeaviside*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter[eN]);
225 epsDirac = epsFactDirac* (useMetrics*h_phi+(1.0-useMetrics)*elementDiameter[eN]);
226 epsDiffusion = epsFactDiffusion*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter[eN]);
238 for(
int i=0;i<nDOF_test_element;i++)
242 register int i_nSpace=i*nSpace;
243 elementResidual_u[i] +=
244 ck.Reaction_weak(
r,u_test_dV[i]) +
245 ck.NumericalDiffusion(epsDiffusion,grad_u,&u_grad_test_dV[i_nSpace]);
248 elementInterface_lumpedMassMatrix[i] += dr*u_test_dV[i];
259 for (
int I=0;I<nSpace;I++)
260 norm += grad_u[I]*grad_u[I];
262 for(
int I=0;I<nSpace;I++)
263 q_n[eN_k_nSpace+I] = grad_u[I]/norm;
269 xt::pyarray<double>& mesh_trial_ref = args.
array<
double>(
"mesh_trial_ref");
270 xt::pyarray<double>& mesh_grad_trial_ref = args.
array<
double>(
"mesh_grad_trial_ref");
271 xt::pyarray<double>& mesh_dof = args.
array<
double>(
"mesh_dof");
272 xt::pyarray<int>& mesh_l2g = args.
array<
int>(
"mesh_l2g");
273 xt::pyarray<double>& dV_ref = args.
array<
double>(
"dV_ref");
274 xt::pyarray<double>& u_trial_ref = args.
array<
double>(
"u_trial_ref");
275 xt::pyarray<double>& u_grad_trial_ref = args.
array<
double>(
"u_grad_trial_ref");
276 xt::pyarray<double>& u_test_ref = args.
array<
double>(
"u_test_ref");
277 xt::pyarray<double>& u_grad_test_ref = args.
array<
double>(
"u_grad_test_ref");
278 xt::pyarray<double>& mesh_trial_trace_ref = args.
array<
double>(
"mesh_trial_trace_ref");
279 xt::pyarray<double>& mesh_grad_trial_trace_ref = args.
array<
double>(
"mesh_grad_trial_trace_ref");
280 xt::pyarray<double>& dS_ref = args.
array<
double>(
"dS_ref");
281 xt::pyarray<double>& u_trial_trace_ref = args.
array<
double>(
"u_trial_trace_ref");
282 xt::pyarray<double>& u_grad_trial_trace_ref = args.
array<
double>(
"u_grad_trial_trace_ref");
283 xt::pyarray<double>& u_test_trace_ref = args.
array<
double>(
"u_test_trace_ref");
284 xt::pyarray<double>& u_grad_test_trace_ref = args.
array<
double>(
"u_grad_test_trace_ref");
285 xt::pyarray<double>& normal_ref = args.
array<
double>(
"normal_ref");
286 xt::pyarray<double>& boundaryJac_ref = args.
array<
double>(
"boundaryJac_ref");
287 int nElements_global = args.
scalar<
int>(
"nElements_global");
288 double useMetrics = args.
scalar<
double>(
"useMetrics");
289 double epsFactHeaviside = args.
scalar<
double>(
"epsFactHeaviside");
290 double epsFactDirac = args.
scalar<
double>(
"epsFactDirac");
291 double epsFactDiffusion = args.
scalar<
double>(
"epsFactDiffusion");
292 xt::pyarray<int>& u_l2g = args.
array<
int>(
"u_l2g");
293 xt::pyarray<double>& elementDiameter = args.
array<
double>(
"elementDiameter");
294 xt::pyarray<double>& nodeDiametersArray = args.
array<
double>(
"nodeDiametersArray");
295 xt::pyarray<double>& u_dof = args.
array<
double>(
"u_dof");
296 xt::pyarray<double>& q_phi = args.
array<
double>(
"q_phi");
297 xt::pyarray<double>& q_normal_phi = args.
array<
double>(
"q_normal_phi");
298 xt::pyarray<double>& ebqe_phi = args.
array<
double>(
"ebqe_phi");
299 xt::pyarray<double>& ebqe_normal_phi = args.
array<
double>(
"ebqe_normal_phi");
300 xt::pyarray<double>& q_H = args.
array<
double>(
"q_H");
301 xt::pyarray<double>& q_u = args.
array<
double>(
"q_u");
302 xt::pyarray<double>& q_n = args.
array<
double>(
"q_n");
303 xt::pyarray<double>& ebqe_u = args.
array<
double>(
"ebqe_u");
304 xt::pyarray<double>& ebqe_n = args.
array<
double>(
"ebqe_n");
305 xt::pyarray<double>& q_r = args.
array<
double>(
"q_r");
306 xt::pyarray<double>& q_vos = args.
array<
double>(
"q_vos");
307 int offset_u = args.
scalar<
int>(
"offset_u");
308 int stride_u = args.
scalar<
int>(
"stride_u");
309 xt::pyarray<double>& globalResidual = args.
array<
double>(
"globalResidual");
310 int nExteriorElementBoundaries_global = args.
scalar<
int>(
"nExteriorElementBoundaries_global");
311 xt::pyarray<int>& exteriorElementBoundariesArray = args.
array<
int>(
"exteriorElementBoundariesArray");
312 xt::pyarray<int>& elementBoundaryElementsArray = args.
array<
int>(
"elementBoundaryElementsArray");
313 xt::pyarray<int>& elementBoundaryLocalElementBoundariesArray = args.
array<
int>(
"elementBoundaryLocalElementBoundariesArray");
314 xt::pyarray<double>& interface_lumpedMassMatrix = args.
array<
double>(
"interface_lumpedMassMatrix");
325 for(
int eN=0;eN<nElements_global;eN++)
328 register double elementResidual_u[nDOF_test_element],
329 element_u[nDOF_trial_element],
330 elementInterface_lumpedMassMatrix[nDOF_test_element];
331 for (
int i=0;i<nDOF_test_element;i++)
333 register int eN_i=eN*nDOF_test_element+i;
334 element_u[i] = u_dof[u_l2g[eN_i]];
337 mesh_grad_trial_ref.data(),
342 u_grad_trial_ref.data(),
344 u_grad_test_ref.data(),
345 mesh_trial_trace_ref.data(),
346 mesh_grad_trial_trace_ref.data(),
348 u_trial_trace_ref.data(),
349 u_grad_trial_trace_ref.data(),
350 u_test_trace_ref.data(),
351 u_grad_test_trace_ref.data(),
353 boundaryJac_ref.data(),
360 elementDiameter.data(),
361 nodeDiametersArray.data(),
366 ebqe_normal_phi.data(),
376 elementInterface_lumpedMassMatrix,
377 nExteriorElementBoundaries_global,
378 exteriorElementBoundariesArray.data(),
379 elementBoundaryElementsArray.data(),
380 elementBoundaryLocalElementBoundariesArray.data(),
386 for(
int i=0;i<nDOF_test_element;i++)
388 register int eN_i=eN*nDOF_test_element+i;
389 int gi = offset_u+stride_u*u_l2g[eN_i];
390 globalResidual[gi] += elementResidual_u[i];
393 interface_lumpedMassMatrix[gi] += elementInterface_lumpedMassMatrix[i];
402 for (
int ebNE = 0; ebNE < nExteriorElementBoundaries_global; ebNE++)
404 register int ebN = exteriorElementBoundariesArray[ebNE],
405 eN = elementBoundaryElementsArray[ebN*2+0],
406 ebN_local = elementBoundaryLocalElementBoundariesArray[ebN*2+0];
409 double element_u[nDOF_trial_element];
410 for (
int i=0;i<nDOF_test_element;i++)
412 register int eN_i=eN*nDOF_test_element+i;
413 element_u[i] = u_dof[u_l2g[eN_i]];
415 for (
int kb=0;kb<nQuadraturePoints_elementBoundary;kb++)
417 register int ebNE_kb = ebNE*nQuadraturePoints_elementBoundary+kb,
418 ebNE_kb_nSpace = ebNE_kb*nSpace,
419 ebN_local_kb = ebN_local*nQuadraturePoints_elementBoundary+kb,
420 ebN_local_kb_nSpace = ebN_local_kb*nSpace;
421 register double u_ext=0.0,
434 jac_ext[nSpace*nSpace],
436 jacInv_ext[nSpace*nSpace],
437 boundaryJac[nSpace*(nSpace-1)],
438 metricTensor[(nSpace-1)*(nSpace-1)],
442 u_grad_trial_trace[nDOF_trial_element*nSpace],
443 normal[nSpace],x_ext,y_ext,z_ext,
444 G[nSpace*nSpace],G_dd_G,tr_G,norm;
448 ck.calculateMapping_elementBoundary(eN,
454 mesh_trial_trace_ref.data(),
455 mesh_grad_trial_trace_ref.data(),
456 boundaryJac_ref.data(),
466 dS = metricTensorDetSqrt*dS_ref[kb];
469 ck.calculateG(jacInv_ext,G,G_dd_G,tr_G);
472 ck.gradTrialFromRef(&u_grad_trial_trace_ref[ebN_local_kb_nSpace*nDOF_trial_element],jacInv_ext,u_grad_trial_trace);
474 ck.valFromElementDOF(element_u,&u_trial_trace_ref[ebN_local_kb*nDOF_test_element],u_ext);
475 ck.gradFromElementDOF(element_u,u_grad_trial_trace,grad_u_ext);
477 ebqe_u[ebNE_kb] = u_ext;
479 for (
int I=0;I<nSpace;I++)
480 norm += grad_u_ext[I]*grad_u_ext[I];
482 for (
int I=0;I<nSpace;I++)
483 ebqe_n[ebNE_kb_nSpace+I] = grad_u_ext[I]/norm;
491 double* mesh_trial_ref,
492 double* mesh_grad_trial_ref,
497 double* u_grad_trial_ref,
499 double* u_grad_test_ref,
501 double* mesh_trial_trace_ref,
502 double* mesh_grad_trial_trace_ref,
504 double* u_trial_trace_ref,
505 double* u_grad_trial_trace_ref,
506 double* u_test_trace_ref,
507 double* u_grad_test_trace_ref,
509 double* boundaryJac_ref,
511 int nElements_global,
513 double epsFactHeaviside,
515 double epsFactDiffusion,
517 double* elementDiameter,
518 double* nodeDiametersArray,
525 double* q_normal_phi,
528 double* elementJacobian_u_u,
532 for (
int i=0;i<nDOF_test_element;i++)
533 for (
int j=0;j<nDOF_trial_element;j++)
535 elementJacobian_u_u[i*nDOF_trial_element+j]=0.0;
537 double epsHeaviside,epsDirac,epsDiffusion;
538 for (
int k=0;k<nQuadraturePoints_element;k++)
540 int eN_k = eN*nQuadraturePoints_element+k,
541 eN_k_nSpace = eN_k*nSpace;
545 register double u=0.0,
550 jacInv[nSpace*nSpace],
551 u_grad_trial[nDOF_trial_element*nSpace],
553 u_test_dV[nDOF_test_element],
554 u_grad_test_dV[nDOF_test_element*nSpace],
556 G[nSpace*nSpace],G_dd_G,tr_G,h_phi;
560 ck.calculateMapping_element(eN,
570 ck.calculateH_element(eN,
577 dV = fabs(jacDet)*dV_ref[k];
578 ck.calculateG(jacInv,G,G_dd_G,tr_G);
591 ck.gradTrialFromRef(&u_grad_trial_ref[k*nDOF_trial_element*nSpace],jacInv,u_grad_trial);
593 ck.valFromElementDOF(element_u,&u_trial_ref[k*nDOF_trial_element],
u);
595 ck.gradFromElementDOF(element_u,u_grad_trial,grad_u);
597 for (
int j=0;j<nDOF_trial_element;j++)
599 u_test_dV[j] = u_test_ref[k*nDOF_trial_element+j]*dV;
600 for (
int I=0;I<nSpace;I++)
602 u_grad_test_dV[j*nSpace+I] = u_grad_trial[j*nSpace+I]*dV;
608 epsHeaviside=epsFactHeaviside*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter[eN]);
609 epsDirac =epsFactDirac* (useMetrics*h_phi+(1.0-useMetrics)*elementDiameter[eN]);
610 epsDiffusion = epsFactDiffusion*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter[eN]);
619 for(
int i=0;i<nDOF_test_element;i++)
623 int i_nSpace=i*nSpace;
624 for(
int j=0;j<nDOF_trial_element;j++)
628 int j_nSpace = j*nSpace;
630 elementJacobian_u_u[i*nDOF_trial_element+j] +=
631 ck.ReactionJacobian_weak(dr,
632 u_trial_ref[k*nDOF_trial_element+j],
634 ck.NumericalDiffusionJacobian(epsDiffusion,
635 &u_grad_trial[j_nSpace],
636 &u_grad_test_dV[i_nSpace]);
643 xt::pyarray<double>& mesh_trial_ref = args.
array<
double>(
"mesh_trial_ref");
644 xt::pyarray<double>& mesh_grad_trial_ref = args.
array<
double>(
"mesh_grad_trial_ref");
645 xt::pyarray<double>& mesh_dof = args.
array<
double>(
"mesh_dof");
646 xt::pyarray<int>& mesh_l2g = args.
array<
int>(
"mesh_l2g");
647 xt::pyarray<double>& dV_ref = args.
array<
double>(
"dV_ref");
648 xt::pyarray<double>& u_trial_ref = args.
array<
double>(
"u_trial_ref");
649 xt::pyarray<double>& u_grad_trial_ref = args.
array<
double>(
"u_grad_trial_ref");
650 xt::pyarray<double>& u_test_ref = args.
array<
double>(
"u_test_ref");
651 xt::pyarray<double>& u_grad_test_ref = args.
array<
double>(
"u_grad_test_ref");
652 xt::pyarray<double>& mesh_trial_trace_ref = args.
array<
double>(
"mesh_trial_trace_ref");
653 xt::pyarray<double>& mesh_grad_trial_trace_ref = args.
array<
double>(
"mesh_grad_trial_trace_ref");
654 xt::pyarray<double>& dS_ref = args.
array<
double>(
"dS_ref");
655 xt::pyarray<double>& u_trial_trace_ref = args.
array<
double>(
"u_trial_trace_ref");
656 xt::pyarray<double>& u_grad_trial_trace_ref = args.
array<
double>(
"u_grad_trial_trace_ref");
657 xt::pyarray<double>& u_test_trace_ref = args.
array<
double>(
"u_test_trace_ref");
658 xt::pyarray<double>& u_grad_test_trace_ref = args.
array<
double>(
"u_grad_test_trace_ref");
659 xt::pyarray<double>& normal_ref = args.
array<
double>(
"normal_ref");
660 xt::pyarray<double>& boundaryJac_ref = args.
array<
double>(
"boundaryJac_ref");
661 int nElements_global = args.
scalar<
int>(
"nElements_global");
662 double useMetrics = args.
scalar<
double>(
"useMetrics");
663 double epsFactHeaviside = args.
scalar<
double>(
"epsFactHeaviside");
664 double epsFactDirac = args.
scalar<
double>(
"epsFactDirac");
665 double epsFactDiffusion = args.
scalar<
double>(
"epsFactDiffusion");
666 xt::pyarray<int>& u_l2g = args.
array<
int>(
"u_l2g");
667 xt::pyarray<double>& elementDiameter = args.
array<
double>(
"elementDiameter");
668 xt::pyarray<double>& nodeDiametersArray = args.
array<
double>(
"nodeDiametersArray");
669 xt::pyarray<double>& u_dof = args.
array<
double>(
"u_dof");
670 xt::pyarray<double>& q_phi = args.
array<
double>(
"q_phi");
671 xt::pyarray<double>& q_normal_phi = args.
array<
double>(
"q_normal_phi");
672 xt::pyarray<double>& q_H = args.
array<
double>(
"q_H");
673 xt::pyarray<double>& q_vos = args.
array<
double>(
"q_vos");
674 xt::pyarray<int>& csrRowIndeces_u_u = args.
array<
int>(
"csrRowIndeces_u_u");
675 xt::pyarray<int>& csrColumnOffsets_u_u = args.
array<
int>(
"csrColumnOffsets_u_u");
676 xt::pyarray<double>& globalJacobian = args.
array<
double>(
"globalJacobian");
677 int numDOFs = args.
scalar<
int>(
"numDOFs");
678 xt::pyarray<int>& csrRowIndeces_DofLoops = args.
array<
int>(
"csrRowIndeces_DofLoops");
679 xt::pyarray<int>& csrColumnOffsets_DofLoops = args.
array<
int>(
"csrColumnOffsets_DofLoops");
680 xt::pyarray<double>& stiffness_matrix = args.
array<
double>(
"stiffness_matrix");
681 xt::pyarray<double>& interface_lumpedMassMatrix = args.
array<
double>(
"interface_lumpedMassMatrix");
685 for (
int i=0; i<numDOFs; i++)
686 for (
int offset=csrRowIndeces_DofLoops[i];
687 offset<csrRowIndeces_DofLoops[i+1]; offset++)
689 int j = csrColumnOffsets_DofLoops[offset];
691 globalJacobian[ij] = interface_lumpedMassMatrix[i] + stiffness_matrix[ij];
693 globalJacobian[ij] = stiffness_matrix[ij];
700 for(
int eN=0;eN<nElements_global;eN++)
702 register double elementJacobian_u_u[nDOF_test_element*nDOF_trial_element],element_u[nDOF_trial_element];
703 for (
int j=0;j<nDOF_trial_element;j++)
705 register int eN_j = eN*nDOF_trial_element+j;
706 element_u[j] = u_dof[u_l2g[eN_j]];
709 mesh_grad_trial_ref.data(),
714 u_grad_trial_ref.data(),
716 u_grad_test_ref.data(),
717 mesh_trial_trace_ref.data(),
718 mesh_grad_trial_trace_ref.data(),
720 u_trial_trace_ref.data(),
721 u_grad_trial_trace_ref.data(),
722 u_test_trace_ref.data(),
723 u_grad_test_trace_ref.data(),
725 boundaryJac_ref.data(),
732 elementDiameter.data(),
733 nodeDiametersArray.data(),
745 for (
int i=0;i<nDOF_test_element;i++)
747 int eN_i = eN*nDOF_test_element+i;
748 for (
int j=0;j<nDOF_trial_element;j++)
750 int eN_i_j = eN_i*nDOF_trial_element+j;
752 globalJacobian[csrRowIndeces_u_u[eN_i] + csrColumnOffsets_u_u[eN_i_j]] += elementJacobian_u_u[i*nDOF_trial_element+j];
761 xt::pyarray<double>& mesh_trial_ref = args.
array<
double>(
"mesh_trial_ref");
762 xt::pyarray<double>& mesh_grad_trial_ref = args.
array<
double>(
"mesh_grad_trial_ref");
763 xt::pyarray<double>& mesh_dof = args.
array<
double>(
"mesh_dof");
764 xt::pyarray<int>& mesh_l2g = args.
array<
int>(
"mesh_l2g");
765 xt::pyarray<double>& dV_ref = args.
array<
double>(
"dV_ref");
766 xt::pyarray<double>& u_trial_ref = args.
array<
double>(
"u_trial_ref");
767 xt::pyarray<double>& u_grad_trial_ref = args.
array<
double>(
"u_grad_trial_ref");
768 xt::pyarray<double>& u_test_ref = args.
array<
double>(
"u_test_ref");
769 xt::pyarray<double>& u_grad_test_ref = args.
array<
double>(
"u_grad_test_ref");
770 xt::pyarray<double>& mesh_trial_trace_ref = args.
array<
double>(
"mesh_trial_trace_ref");
771 xt::pyarray<double>& mesh_grad_trial_trace_ref = args.
array<
double>(
"mesh_grad_trial_trace_ref");
772 xt::pyarray<double>& dS_ref = args.
array<
double>(
"dS_ref");
773 xt::pyarray<double>& u_trial_trace_ref = args.
array<
double>(
"u_trial_trace_ref");
774 xt::pyarray<double>& u_grad_trial_trace_ref = args.
array<
double>(
"u_grad_trial_trace_ref");
775 xt::pyarray<double>& u_test_trace_ref = args.
array<
double>(
"u_test_trace_ref");
776 xt::pyarray<double>& u_grad_test_trace_ref = args.
array<
double>(
"u_grad_test_trace_ref");
777 xt::pyarray<double>& normal_ref = args.
array<
double>(
"normal_ref");
778 xt::pyarray<double>& boundaryJac_ref = args.
array<
double>(
"boundaryJac_ref");
779 int nElements_global = args.
scalar<
int>(
"nElements_global");
780 double useMetrics = args.
scalar<
double>(
"useMetrics");
781 double epsFactHeaviside = args.
scalar<
double>(
"epsFactHeaviside");
782 double epsFactDirac = args.
scalar<
double>(
"epsFactDirac");
783 double epsFactDiffusion = args.
scalar<
double>(
"epsFactDiffusion");
784 xt::pyarray<int>& u_l2g = args.
array<
int>(
"u_l2g");
785 xt::pyarray<double>& elementDiameter = args.
array<
double>(
"elementDiameter");
786 xt::pyarray<double>& nodeDiametersArray = args.
array<
double>(
"nodeDiametersArray");
787 xt::pyarray<double>& u_dof = args.
array<
double>(
"u_dof");
788 xt::pyarray<double>& q_phi = args.
array<
double>(
"q_phi");
789 xt::pyarray<double>& q_normal_phi = args.
array<
double>(
"q_normal_phi");
790 xt::pyarray<double>& ebqe_phi = args.
array<
double>(
"ebqe_phi");
791 xt::pyarray<double>& ebqe_normal_phi = args.
array<
double>(
"ebqe_normal_phi");
792 xt::pyarray<double>& q_H = args.
array<
double>(
"q_H");
793 xt::pyarray<double>& q_u = args.
array<
double>(
"q_u");
794 xt::pyarray<double>& q_n = args.
array<
double>(
"q_n");
795 xt::pyarray<double>& ebqe_u = args.
array<
double>(
"ebqe_u");
796 xt::pyarray<double>& ebqe_n = args.
array<
double>(
"ebqe_n");
797 xt::pyarray<double>& q_r = args.
array<
double>(
"q_r");
798 xt::pyarray<double>& q_vos = args.
array<
double>(
"q_vos");
799 int offset_u = args.
scalar<
int>(
"offset_u");
800 int stride_u = args.
scalar<
int>(
"stride_u");
801 xt::pyarray<double>& globalResidual = args.
array<
double>(
"globalResidual");
802 int nExteriorElementBoundaries_global = args.
scalar<
int>(
"nExteriorElementBoundaries_global");
803 xt::pyarray<int>& exteriorElementBoundariesArray = args.
array<
int>(
"exteriorElementBoundariesArray");
804 xt::pyarray<int>& elementBoundaryElementsArray = args.
array<
int>(
"elementBoundaryElementsArray");
805 xt::pyarray<int>& elementBoundaryLocalElementBoundariesArray = args.
array<
int>(
"elementBoundaryLocalElementBoundariesArray");
806 int maxIts = args.
scalar<
int>(
"maxIts");
807 double atol = args.
scalar<
double>(
"atol");
818 for(
int eN=0;eN<nElements_global;eN++)
821 register double element_u[nDOF_test_element],
822 element_du[nDOF_test_element],
823 elementResidual_u[nDOF_test_element],
824 dummy[nDOF_test_element],
825 elementJacobian_u_u[nDOF_test_element*nDOF_trial_element],scale=1.0;
826 register PROTEUS_LAPACK_INTEGER elementPivots[nDOF_test_element],
827 elementColPivots[nDOF_test_element];
829 for (
int i=0;i<nDOF_test_element;i++)
834 mesh_grad_trial_ref.data(),
839 u_grad_trial_ref.data(),
841 u_grad_test_ref.data(),
842 mesh_trial_trace_ref.data(),
843 mesh_grad_trial_trace_ref.data(),
845 u_trial_trace_ref.data(),
846 u_grad_trial_trace_ref.data(),
847 u_test_trace_ref.data(),
848 u_grad_test_trace_ref.data(),
850 boundaryJac_ref.data(),
857 elementDiameter.data(),
858 nodeDiametersArray.data(),
863 ebqe_normal_phi.data(),
874 nExteriorElementBoundaries_global,
875 exteriorElementBoundariesArray.data(),
876 elementBoundaryElementsArray.data(),
877 elementBoundaryLocalElementBoundariesArray.data(),
882 for (
int i=0;i<nDOF_test_element;i++)
884 resNorm += elementResidual_u[i];
886 resNorm = fabs(resNorm);
891 while (resNorm >= atol && its < maxIts)
895 mesh_grad_trial_ref.data(),
900 u_grad_trial_ref.data(),
902 u_grad_test_ref.data(),
903 mesh_trial_trace_ref.data(),
904 mesh_grad_trial_trace_ref.data(),
906 u_trial_trace_ref.data(),
907 u_grad_trial_trace_ref.data(),
908 u_test_trace_ref.data(),
909 u_grad_test_trace_ref.data(),
911 boundaryJac_ref.data(),
918 elementDiameter.data(),
919 nodeDiametersArray.data(),
928 for (
int i=0;i<nDOF_test_element;i++)
930 element_du[i] = -elementResidual_u[i];
931 elementPivots[i] = ((PROTEUS_LAPACK_INTEGER)0);
932 elementColPivots[i]=((PROTEUS_LAPACK_INTEGER)0);
941 PROTEUS_LAPACK_INTEGER La_N=((PROTEUS_LAPACK_INTEGER)nDOF_test_element),
957 double resNormNew = resNorm,lambda=1.0;
959 while (resNormNew > 0.99*resNorm && lsIts < 100)
962 for (
int i=0;i<nDOF_test_element;i++)
964 element_u[i] += lambda*element_du[i];
969 mesh_grad_trial_ref.data(),
974 u_grad_trial_ref.data(),
976 u_grad_test_ref.data(),
977 mesh_trial_trace_ref.data(),
978 mesh_grad_trial_trace_ref.data(),
980 u_trial_trace_ref.data(),
981 u_grad_trial_trace_ref.data(),
982 u_test_trace_ref.data(),
983 u_grad_test_trace_ref.data(),
985 boundaryJac_ref.data(),
992 elementDiameter.data(),
993 nodeDiametersArray.data(),
998 ebqe_normal_phi.data(),
1009 nExteriorElementBoundaries_global,
1010 exteriorElementBoundariesArray.data(),
1011 elementBoundaryElementsArray.data(),
1012 elementBoundaryLocalElementBoundariesArray.data(),
1018 for (
int i=0;i<nDOF_test_element;i++)
1020 resNormNew += elementResidual_u[i];
1021 std::cout<<
"element_u["<<i<<
"] "<<element_u[i]<<std::endl;
1022 std::cout<<
"elementResidual_u["<<i<<
"] "<<elementResidual_u[i]<<std::endl;
1024 resNormNew = fabs(resNormNew);
1026 resNorm = resNormNew;
1027 std::cout<<
"INFO "<<INFO<<std::endl;
1028 std::cout<<
"resNorm["<<its<<
"] "<<resNorm<<std::endl;
1034 xt::pyarray<double>& mesh_trial_ref = args.
array<
double>(
"mesh_trial_ref");
1035 xt::pyarray<double>& mesh_grad_trial_ref = args.
array<
double>(
"mesh_grad_trial_ref");
1036 xt::pyarray<double>& mesh_dof = args.
array<
double>(
"mesh_dof");
1037 xt::pyarray<int>& mesh_l2g = args.
array<
int>(
"mesh_l2g");
1038 xt::pyarray<double>& dV_ref = args.
array<
double>(
"dV_ref");
1039 xt::pyarray<double>& u_trial_ref = args.
array<
double>(
"u_trial_ref");
1040 xt::pyarray<double>& u_grad_trial_ref = args.
array<
double>(
"u_grad_trial_ref");
1041 xt::pyarray<double>& u_test_ref = args.
array<
double>(
"u_test_ref");
1042 xt::pyarray<double>& u_grad_test_ref = args.
array<
double>(
"u_grad_test_ref");
1043 xt::pyarray<double>& mesh_trial_trace_ref = args.
array<
double>(
"mesh_trial_trace_ref");
1044 xt::pyarray<double>& mesh_grad_trial_trace_ref = args.
array<
double>(
"mesh_grad_trial_trace_ref");
1045 xt::pyarray<double>& dS_ref = args.
array<
double>(
"dS_ref");
1046 xt::pyarray<double>& u_trial_trace_ref = args.
array<
double>(
"u_trial_trace_ref");
1047 xt::pyarray<double>& u_grad_trial_trace_ref = args.
array<
double>(
"u_grad_trial_trace_ref");
1048 xt::pyarray<double>& u_test_trace_ref = args.
array<
double>(
"u_test_trace_ref");
1049 xt::pyarray<double>& u_grad_test_trace_ref = args.
array<
double>(
"u_grad_test_trace_ref");
1050 xt::pyarray<double>& normal_ref = args.
array<
double>(
"normal_ref");
1051 xt::pyarray<double>& boundaryJac_ref = args.
array<
double>(
"boundaryJac_ref");
1052 int nElements_global = args.
scalar<
int>(
"nElements_global");
1053 double useMetrics = args.
scalar<
double>(
"useMetrics");
1054 double epsFactHeaviside = args.
scalar<
double>(
"epsFactHeaviside");
1055 double epsFactDirac = args.
scalar<
double>(
"epsFactDirac");
1056 double epsFactDiffusion = args.
scalar<
double>(
"epsFactDiffusion");
1057 xt::pyarray<int>& u_l2g = args.
array<
int>(
"u_l2g");
1058 xt::pyarray<double>& elementDiameter = args.
array<
double>(
"elementDiameter");
1059 xt::pyarray<double>& nodeDiametersArray = args.
array<
double>(
"nodeDiametersArray");
1060 xt::pyarray<double>& u_dof = args.
array<
double>(
"u_dof");
1061 xt::pyarray<double>& q_phi = args.
array<
double>(
"q_phi");
1062 xt::pyarray<double>& q_normal_phi = args.
array<
double>(
"q_normal_phi");
1063 xt::pyarray<double>& ebqe_phi = args.
array<
double>(
"ebqe_phi");
1064 xt::pyarray<double>& ebqe_normal_phi = args.
array<
double>(
"ebqe_normal_phi");
1065 xt::pyarray<double>& q_H = args.
array<
double>(
"q_H");
1066 xt::pyarray<double>& q_u = args.
array<
double>(
"q_u");
1067 xt::pyarray<double>& q_n = args.
array<
double>(
"q_n");
1068 xt::pyarray<double>& ebqe_u = args.
array<
double>(
"ebqe_u");
1069 xt::pyarray<double>& ebqe_n = args.
array<
double>(
"ebqe_n");
1070 xt::pyarray<double>& q_r = args.
array<
double>(
"q_r");
1071 xt::pyarray<double>& q_vos = args.
array<
double>(
"q_vos");
1072 int offset_u = args.
scalar<
int>(
"offset_u");
1073 int stride_u = args.
scalar<
int>(
"stride_u");
1074 xt::pyarray<double>& globalResidual = args.
array<
double>(
"globalResidual");
1075 int nExteriorElementBoundaries_global = args.
scalar<
int>(
"nExteriorElementBoundaries_global");
1076 xt::pyarray<int>& exteriorElementBoundariesArray = args.
array<
int>(
"exteriorElementBoundariesArray");
1077 xt::pyarray<int>& elementBoundaryElementsArray = args.
array<
int>(
"elementBoundaryElementsArray");
1078 xt::pyarray<int>& elementBoundaryLocalElementBoundariesArray = args.
array<
int>(
"elementBoundaryLocalElementBoundariesArray");
1079 int maxIts = args.
scalar<
int>(
"maxIts");
1080 double atol = args.
scalar<
double>(
"atol");
1081 for(
int eN=0;eN<nElements_global;eN++)
1084 register double element_u[nDOF_test_element],elementConstant_u,
1085 elementResidual_u[nDOF_test_element],
1086 dummy[nDOF_test_element],
1087 elementConstantResidual,
1088 elementJacobian_u_u[nDOF_test_element*nDOF_trial_element],elementConstantJacobian,resNorm;
1089 elementConstant_u=0.0;
1090 for (
int i=0;i<nDOF_test_element;i++)
1092 element_u[i]=elementConstant_u;
1095 mesh_grad_trial_ref.data(),
1100 u_grad_trial_ref.data(),
1102 u_grad_test_ref.data(),
1103 mesh_trial_trace_ref.data(),
1104 mesh_grad_trial_trace_ref.data(),
1106 u_trial_trace_ref.data(),
1107 u_grad_trial_trace_ref.data(),
1108 u_test_trace_ref.data(),
1109 u_grad_test_trace_ref.data(),
1111 boundaryJac_ref.data(),
1118 elementDiameter.data(),
1119 nodeDiametersArray.data(),
1122 q_normal_phi.data(),
1124 ebqe_normal_phi.data(),
1135 nExteriorElementBoundaries_global,
1136 exteriorElementBoundariesArray.data(),
1137 elementBoundaryElementsArray.data(),
1138 elementBoundaryLocalElementBoundariesArray.data(),
1142 elementConstantResidual=0.0;
1143 for (
int i=0;i<nDOF_test_element;i++)
1145 elementConstantResidual += elementResidual_u[i];
1147 resNorm = fabs(elementConstantResidual);
1152 while (resNorm >= atol && its < maxIts)
1156 mesh_grad_trial_ref.data(),
1161 u_grad_trial_ref.data(),
1163 u_grad_test_ref.data(),
1164 mesh_trial_trace_ref.data(),
1165 mesh_grad_trial_trace_ref.data(),
1167 u_trial_trace_ref.data(),
1168 u_grad_trial_trace_ref.data(),
1169 u_test_trace_ref.data(),
1170 u_grad_test_trace_ref.data(),
1172 boundaryJac_ref.data(),
1179 elementDiameter.data(),
1180 nodeDiametersArray.data(),
1183 q_normal_phi.data(),
1186 elementJacobian_u_u,
1189 elementConstantJacobian=0.0;
1190 for (
int i=0;i<nDOF_test_element;i++)
1192 for (
int j=0;j<nDOF_test_element;j++)
1194 elementConstantJacobian += elementJacobian_u_u[i*nDOF_trial_element+j];
1197 std::cout<<
"elementConstantJacobian "<<elementConstantJacobian<<std::endl;
1199 elementConstant_u -= elementConstantResidual/(elementConstantJacobian+1.0e-8);
1200 for (
int i=0;i<nDOF_test_element;i++)
1202 element_u[i] = elementConstant_u;
1206 mesh_grad_trial_ref.data(),
1211 u_grad_trial_ref.data(),
1213 u_grad_test_ref.data(),
1214 mesh_trial_trace_ref.data(),
1215 mesh_grad_trial_trace_ref.data(),
1217 u_trial_trace_ref.data(),
1218 u_grad_trial_trace_ref.data(),
1219 u_test_trace_ref.data(),
1220 u_grad_test_trace_ref.data(),
1222 boundaryJac_ref.data(),
1229 elementDiameter.data(),
1230 nodeDiametersArray.data(),
1233 q_normal_phi.data(),
1235 ebqe_normal_phi.data(),
1246 nExteriorElementBoundaries_global,
1247 exteriorElementBoundariesArray.data(),
1248 elementBoundaryElementsArray.data(),
1249 elementBoundaryLocalElementBoundariesArray.data(),
1253 elementConstantResidual=0.0;
1254 for (
int i=0;i<nDOF_test_element;i++)
1256 elementConstantResidual += elementResidual_u[i];
1258 resNorm = fabs(elementConstantResidual);
1259 std::cout<<
"resNorm["<<its<<
"] "<<resNorm<<std::endl;
1266 xt::pyarray<double>& mesh_trial_ref = args.
array<
double>(
"mesh_trial_ref");
1267 xt::pyarray<double>& mesh_grad_trial_ref = args.
array<
double>(
"mesh_grad_trial_ref");
1268 xt::pyarray<double>& mesh_dof = args.
array<
double>(
"mesh_dof");
1269 xt::pyarray<int>& mesh_l2g = args.
array<
int>(
"mesh_l2g");
1270 xt::pyarray<double>& dV_ref = args.
array<
double>(
"dV_ref");
1271 xt::pyarray<double>& u_trial_ref = args.
array<
double>(
"u_trial_ref");
1272 xt::pyarray<double>& u_grad_trial_ref = args.
array<
double>(
"u_grad_trial_ref");
1273 xt::pyarray<double>& u_test_ref = args.
array<
double>(
"u_test_ref");
1274 xt::pyarray<double>& u_grad_test_ref = args.
array<
double>(
"u_grad_test_ref");
1275 xt::pyarray<double>& mesh_trial_trace_ref = args.
array<
double>(
"mesh_trial_trace_ref");
1276 xt::pyarray<double>& mesh_grad_trial_trace_ref = args.
array<
double>(
"mesh_grad_trial_trace_ref");
1277 xt::pyarray<double>& dS_ref = args.
array<
double>(
"dS_ref");
1278 xt::pyarray<double>& u_trial_trace_ref = args.
array<
double>(
"u_trial_trace_ref");
1279 xt::pyarray<double>& u_grad_trial_trace_ref = args.
array<
double>(
"u_grad_trial_trace_ref");
1280 xt::pyarray<double>& u_test_trace_ref = args.
array<
double>(
"u_test_trace_ref");
1281 xt::pyarray<double>& u_grad_test_trace_ref = args.
array<
double>(
"u_grad_test_trace_ref");
1282 xt::pyarray<double>& normal_ref = args.
array<
double>(
"normal_ref");
1283 xt::pyarray<double>& boundaryJac_ref = args.
array<
double>(
"boundaryJac_ref");
1284 int nElements_owned = args.
scalar<
int>(
"nElements_owned");
1285 double useMetrics = args.
scalar<
double>(
"useMetrics");
1286 double epsFactHeaviside = args.
scalar<
double>(
"epsFactHeaviside");
1287 double epsFactDirac = args.
scalar<
double>(
"epsFactDirac");
1288 double epsFactDiffusion = args.
scalar<
double>(
"epsFactDiffusion");
1289 xt::pyarray<int>& u_l2g = args.
array<
int>(
"u_l2g");
1290 xt::pyarray<double>& elementDiameter = args.
array<
double>(
"elementDiameter");
1291 xt::pyarray<double>& nodeDiametersArray = args.
array<
double>(
"nodeDiametersArray");
1292 xt::pyarray<double>& u_dof = args.
array<
double>(
"u_dof");
1293 xt::pyarray<double>& q_phi = args.
array<
double>(
"q_phi");
1294 xt::pyarray<double>& q_normal_phi = args.
array<
double>(
"q_normal_phi");
1295 xt::pyarray<double>& ebqe_phi = args.
array<
double>(
"ebqe_phi");
1296 xt::pyarray<double>& ebqe_normal_phi = args.
array<
double>(
"ebqe_normal_phi");
1297 xt::pyarray<double>& q_H = args.
array<
double>(
"q_H");
1298 xt::pyarray<double>& q_u = args.
array<
double>(
"q_u");
1299 xt::pyarray<double>& q_n = args.
array<
double>(
"q_n");
1300 xt::pyarray<double>& ebqe_u = args.
array<
double>(
"ebqe_u");
1301 xt::pyarray<double>& ebqe_n = args.
array<
double>(
"ebqe_n");
1302 xt::pyarray<double>& q_r = args.
array<
double>(
"q_r");
1303 xt::pyarray<double>& q_vos = args.
array<
double>(
"q_vos");
1304 int offset_u = args.
scalar<
int>(
"offset_u");
1305 int stride_u = args.
scalar<
int>(
"stride_u");
1306 xt::pyarray<double>& globalResidual = args.
array<
double>(
"globalResidual");
1307 int nExteriorElementBoundaries_global = args.
scalar<
int>(
"nExteriorElementBoundaries_global");
1308 xt::pyarray<int>& exteriorElementBoundariesArray = args.
array<
int>(
"exteriorElementBoundariesArray");
1309 xt::pyarray<int>& elementBoundaryElementsArray = args.
array<
int>(
"elementBoundaryElementsArray");
1310 xt::pyarray<int>& elementBoundaryLocalElementBoundariesArray = args.
array<
int>(
"elementBoundaryLocalElementBoundariesArray");
1311 int maxIts = args.
scalar<
int>(
"maxIts");
1312 double atol = args.
scalar<
double>(
"atol");
1313 double constant_u = args.
scalar<
double>(
"constant_u");
1314 register double element_u[nDOF_test_element],
1315 elementResidual_u[nDOF_test_element],
1316 dummy[nDOF_test_element],
1317 elementJacobian_u_u[nDOF_test_element*nDOF_trial_element];
1318 double constantResidual = 0.0;
1319 double constantJacobian = 0.0;
1320 for (
int i=0;i<nDOF_trial_element;i++)
1322 element_u[i]=constant_u;
1325 for(
int eN=0;eN<nElements_owned;eN++)
1328 mesh_grad_trial_ref.data(),
1333 u_grad_trial_ref.data(),
1335 u_grad_test_ref.data(),
1336 mesh_trial_trace_ref.data(),
1337 mesh_grad_trial_trace_ref.data(),
1339 u_trial_trace_ref.data(),
1340 u_grad_trial_trace_ref.data(),
1341 u_test_trace_ref.data(),
1342 u_grad_test_trace_ref.data(),
1344 boundaryJac_ref.data(),
1351 elementDiameter.data(),
1352 nodeDiametersArray.data(),
1355 q_normal_phi.data(),
1357 ebqe_normal_phi.data(),
1368 nExteriorElementBoundaries_global,
1369 exteriorElementBoundariesArray.data(),
1370 elementBoundaryElementsArray.data(),
1371 elementBoundaryLocalElementBoundariesArray.data(),
1375 for (
int i=0;i<nDOF_test_element;i++)
1377 constantResidual += elementResidual_u[i];
1380 mesh_grad_trial_ref.data(),
1385 u_grad_trial_ref.data(),
1387 u_grad_test_ref.data(),
1388 mesh_trial_trace_ref.data(),
1389 mesh_grad_trial_trace_ref.data(),
1391 u_trial_trace_ref.data(),
1392 u_grad_trial_trace_ref.data(),
1393 u_test_trace_ref.data(),
1394 u_grad_test_trace_ref.data(),
1396 boundaryJac_ref.data(),
1403 elementDiameter.data(),
1404 nodeDiametersArray.data(),
1407 q_normal_phi.data(),
1410 elementJacobian_u_u,
1413 for (
int i=0;i<nDOF_test_element;i++)
1415 for (
int j=0;j<nDOF_test_element;j++)
1417 constantJacobian += elementJacobian_u_u[i*nDOF_trial_element+j];
1421 return std::make_pair(constantResidual, constantJacobian);
1426 xt::pyarray<double>& mesh_trial_ref = args.
array<
double>(
"mesh_trial_ref");
1427 xt::pyarray<double>& mesh_grad_trial_ref = args.
array<
double>(
"mesh_grad_trial_ref");
1428 xt::pyarray<double>& mesh_dof = args.
array<
double>(
"mesh_dof");
1429 xt::pyarray<int>& mesh_l2g = args.
array<
int>(
"mesh_l2g");
1430 xt::pyarray<double>& dV_ref = args.
array<
double>(
"dV_ref");
1431 xt::pyarray<double>& u_trial_ref = args.
array<
double>(
"u_trial_ref");
1432 xt::pyarray<double>& u_grad_trial_ref = args.
array<
double>(
"u_grad_trial_ref");
1433 xt::pyarray<double>& u_test_ref = args.
array<
double>(
"u_test_ref");
1434 xt::pyarray<double>& u_grad_test_ref = args.
array<
double>(
"u_grad_test_ref");
1435 xt::pyarray<double>& mesh_trial_trace_ref = args.
array<
double>(
"mesh_trial_trace_ref");
1436 xt::pyarray<double>& mesh_grad_trial_trace_ref = args.
array<
double>(
"mesh_grad_trial_trace_ref");
1437 xt::pyarray<double>& dS_ref = args.
array<
double>(
"dS_ref");
1438 xt::pyarray<double>& u_trial_trace_ref = args.
array<
double>(
"u_trial_trace_ref");
1439 xt::pyarray<double>& u_grad_trial_trace_ref = args.
array<
double>(
"u_grad_trial_trace_ref");
1440 xt::pyarray<double>& u_test_trace_ref = args.
array<
double>(
"u_test_trace_ref");
1441 xt::pyarray<double>& u_grad_test_trace_ref = args.
array<
double>(
"u_grad_test_trace_ref");
1442 xt::pyarray<double>& normal_ref = args.
array<
double>(
"normal_ref");
1443 xt::pyarray<double>& boundaryJac_ref = args.
array<
double>(
"boundaryJac_ref");
1444 int nElements_owned = args.
scalar<
int>(
"nElements_owned");
1445 double useMetrics = args.
scalar<
double>(
"useMetrics");
1446 double epsFactHeaviside = args.
scalar<
double>(
"epsFactHeaviside");
1447 double epsFactDirac = args.
scalar<
double>(
"epsFactDirac");
1448 double epsFactDiffusion = args.
scalar<
double>(
"epsFactDiffusion");
1449 xt::pyarray<int>& u_l2g = args.
array<
int>(
"u_l2g");
1450 xt::pyarray<double>& elementDiameter = args.
array<
double>(
"elementDiameter");
1451 xt::pyarray<double>& nodeDiametersArray = args.
array<
double>(
"nodeDiametersArray");
1452 xt::pyarray<double>& u_dof = args.
array<
double>(
"u_dof");
1453 xt::pyarray<double>& q_phi = args.
array<
double>(
"q_phi");
1454 xt::pyarray<double>& q_normal_phi = args.
array<
double>(
"q_normal_phi");
1455 xt::pyarray<double>& ebqe_phi = args.
array<
double>(
"ebqe_phi");
1456 xt::pyarray<double>& ebqe_normal_phi = args.
array<
double>(
"ebqe_normal_phi");
1457 xt::pyarray<double>& q_H = args.
array<
double>(
"q_H");
1458 xt::pyarray<double>& q_u = args.
array<
double>(
"q_u");
1459 xt::pyarray<double>& q_n = args.
array<
double>(
"q_n");
1460 xt::pyarray<double>& ebqe_u = args.
array<
double>(
"ebqe_u");
1461 xt::pyarray<double>& ebqe_n = args.
array<
double>(
"ebqe_n");
1462 xt::pyarray<double>& q_r = args.
array<
double>(
"q_r");
1463 xt::pyarray<double>& q_vos = args.
array<
double>(
"q_vos");
1464 int offset_u = args.
scalar<
int>(
"offset_u");
1465 int stride_u = args.
scalar<
int>(
"stride_u");
1466 xt::pyarray<double>& globalResidual = args.
array<
double>(
"globalResidual");
1467 int nExteriorElementBoundaries_global = args.
scalar<
int>(
"nExteriorElementBoundaries_global");
1468 xt::pyarray<int>& exteriorElementBoundariesArray = args.
array<
int>(
"exteriorElementBoundariesArray");
1469 xt::pyarray<int>& elementBoundaryElementsArray = args.
array<
int>(
"elementBoundaryElementsArray");
1470 xt::pyarray<int>& elementBoundaryLocalElementBoundariesArray = args.
array<
int>(
"elementBoundaryLocalElementBoundariesArray");
1471 double globalMass = 0.0;
1472 for(
int eN=0;eN<nElements_owned;eN++)
1474 double epsHeaviside;
1476 for (
int k=0;k<nQuadraturePoints_element;k++)
1479 register int eN_k = eN*nQuadraturePoints_element+k,
1480 eN_k_nSpace = eN_k*nSpace;
1483 double jac[nSpace*nSpace],
1485 jacInv[nSpace*nSpace],
1490 G[nSpace*nSpace],G_dd_G,tr_G,h_phi;
1494 ck.calculateMapping_element(eN,
1498 mesh_trial_ref.data(),
1499 mesh_grad_trial_ref.data(),
1504 ck.calculateH_element(eN,
1506 nodeDiametersArray.data(),
1508 mesh_trial_ref.data(),
1511 dV = fabs(jacDet)*dV_ref[k];
1512 ck.calculateG(jacInv,G,G_dd_G,tr_G);
1521 epsHeaviside=epsFactHeaviside*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter[eN]);
1529 xt::pyarray<double>& mesh_trial_ip = args.
array<
double>(
"mesh_trial_ip");
1530 xt::pyarray<double>& mesh_trial_ref = args.
array<
double>(
"mesh_trial_ref");
1531 xt::pyarray<double>& mesh_grad_trial_ref = args.
array<
double>(
"mesh_grad_trial_ref");
1532 xt::pyarray<double>& mesh_dof = args.
array<
double>(
"mesh_dof");
1533 xt::pyarray<int>& mesh_l2g = args.
array<
int>(
"mesh_l2g");
1534 xt::pyarray<double>& dV_ref = args.
array<
double>(
"dV_ref");
1535 xt::pyarray<double>& u_trial_ref = args.
array<
double>(
"u_trial_ref");
1536 xt::pyarray<double>& u_grad_trial_ref = args.
array<
double>(
"u_grad_trial_ref");
1537 xt::pyarray<double>& u_test_ref = args.
array<
double>(
"u_test_ref");
1538 xt::pyarray<double>& u_grad_test_ref = args.
array<
double>(
"u_grad_test_ref");
1539 xt::pyarray<double>& mesh_trial_trace_ref = args.
array<
double>(
"mesh_trial_trace_ref");
1540 xt::pyarray<double>& mesh_grad_trial_trace_ref = args.
array<
double>(
"mesh_grad_trial_trace_ref");
1541 xt::pyarray<double>& dS_ref = args.
array<
double>(
"dS_ref");
1542 xt::pyarray<double>& u_trial_trace_ref = args.
array<
double>(
"u_trial_trace_ref");
1543 xt::pyarray<double>& u_grad_trial_trace_ref = args.
array<
double>(
"u_grad_trial_trace_ref");
1544 xt::pyarray<double>& u_test_trace_ref = args.
array<
double>(
"u_test_trace_ref");
1545 xt::pyarray<double>& u_grad_test_trace_ref = args.
array<
double>(
"u_grad_test_trace_ref");
1546 xt::pyarray<double>& normal_ref = args.
array<
double>(
"normal_ref");
1547 xt::pyarray<double>& boundaryJac_ref = args.
array<
double>(
"boundaryJac_ref");
1548 int nElements_global = args.
scalar<
int>(
"nElements_global");
1549 double useMetrics = args.
scalar<
double>(
"useMetrics");
1550 double epsFactHeaviside = args.
scalar<
double>(
"epsFactHeaviside");
1551 double epsFactDirac = args.
scalar<
double>(
"epsFactDirac");
1552 double epsFactDiffusion = args.
scalar<
double>(
"epsFactDiffusion");
1553 xt::pyarray<int>& phi_l2g = args.
array<
int>(
"phi_l2g");
1554 xt::pyarray<double>& elementDiameter = args.
array<
double>(
"elementDiameter");
1555 xt::pyarray<double>& nodeDiametersArray = args.
array<
double>(
"nodeDiametersArray");
1556 xt::pyarray<double>& phi_dof = args.
array<
double>(
"phi_dof");
1557 xt::pyarray<double>& q_phi = args.
array<
double>(
"q_phi");
1558 xt::pyarray<double>& q_normal_phi = args.
array<
double>(
"q_normal_phi");
1559 xt::pyarray<double>& ebqe_phi = args.
array<
double>(
"ebqe_phi");
1560 xt::pyarray<double>& ebqe_normal_phi = args.
array<
double>(
"ebqe_normal_phi");
1561 xt::pyarray<double>& q_H = args.
array<
double>(
"q_H");
1562 xt::pyarray<double>& q_u = args.
array<
double>(
"q_u");
1563 xt::pyarray<double>& q_n = args.
array<
double>(
"q_n");
1564 xt::pyarray<double>& ebqe_u = args.
array<
double>(
"ebqe_u");
1565 xt::pyarray<double>& ebqe_n = args.
array<
double>(
"ebqe_n");
1566 xt::pyarray<double>& q_r = args.
array<
double>(
"q_r");
1567 xt::pyarray<double>& q_vos = args.
array<
double>(
"q_vos");
1568 int offset_u = args.
scalar<
int>(
"offset_u");
1569 int stride_u = args.
scalar<
int>(
"stride_u");
1570 xt::pyarray<double>& globalResidual = args.
array<
double>(
"globalResidual");
1571 int nExteriorElementBoundaries_global = args.
scalar<
int>(
"nExteriorElementBoundaries_global");
1572 xt::pyarray<int>& exteriorElementBoundariesArray = args.
array<
int>(
"exteriorElementBoundariesArray");
1573 xt::pyarray<int>& elementBoundaryElementsArray = args.
array<
int>(
"elementBoundaryElementsArray");
1574 xt::pyarray<int>& elementBoundaryLocalElementBoundariesArray = args.
array<
int>(
"elementBoundaryLocalElementBoundariesArray");
1575 xt::pyarray<double>& H_dof = args.
array<
double>(
"H_dof");
1576 for(
int eN=0;eN<nElements_global;eN++)
1578 double epsHeaviside;
1580 for (
int k=0;k<nQuadraturePoints_element;k++)
1583 register int eN_k = eN*nQuadraturePoints_element+k,
1584 eN_k_nSpace = eN_k*nSpace;
1587 register double jac[nSpace*nSpace],
1589 jacInv[nSpace*nSpace],
1594 G[nSpace*nSpace],G_dd_G,tr_G,h_phi;
1598 ck.calculateMapping_element(eN,
1602 mesh_trial_ref.data(),
1603 mesh_grad_trial_ref.data(),
1608 ck.calculateH_element(eN,
1610 nodeDiametersArray.data(),
1612 mesh_trial_ref.data(),
1615 dV = fabs(jacDet)*dV_ref[k];
1616 ck.calculateG(jacInv,G,G_dd_G,tr_G);
1626 epsHeaviside=epsFactHeaviside*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter[eN]);
1629 for (
int i=0;i<nDOF_trial_element;i++)
1631 int eN_i = eN*nDOF_trial_element + i;
1632 register double h_phi=0.0;
1633 ck.calculateH_element(eN,
1635 nodeDiametersArray.data(),
1637 mesh_trial_ip.data(),
1639 epsHeaviside = epsFactHeaviside*h_phi;
1647 xt::pyarray<double>& mesh_trial_ref = args.
array<
double>(
"mesh_trial_ref");
1648 xt::pyarray<double>& mesh_grad_trial_ref = args.
array<
double>(
"mesh_grad_trial_ref");
1649 xt::pyarray<double>& mesh_dof = args.
array<
double>(
"mesh_dof");
1650 xt::pyarray<int>& mesh_l2g = args.
array<
int>(
"mesh_l2g");
1651 xt::pyarray<double>& dV_ref = args.
array<
double>(
"dV_ref");
1652 xt::pyarray<double>& u_grad_trial_ref = args.
array<
double>(
"u_grad_trial_ref");
1653 int nElements_global = args.
scalar<
int>(
"nElements_global");
1654 xt::pyarray<int>& csrRowIndeces_u_u = args.
array<
int>(
"csrRowIndeces_u_u");
1655 xt::pyarray<int>& csrColumnOffsets_u_u = args.
array<
int>(
"csrColumnOffsets_u_u");
1656 xt::pyarray<double>& globalJacobian = args.
array<
double>(
"globalJacobian");
1657 double useMetrics = args.
scalar<
double>(
"useMetrics");
1658 double epsFactDiffusion = args.
scalar<
double>(
"epsFactDiffusion");
1659 xt::pyarray<double>& elementDiameter = args.
array<
double>(
"elementDiameter");
1660 xt::pyarray<double>& nodeDiametersArray = args.
array<
double>(
"nodeDiametersArray");
1664 for(
int eN=0;eN<nElements_global;eN++)
1666 register double elementJacobian_u_u[nDOF_test_element][nDOF_trial_element];
1667 for (
int i=0;i<nDOF_test_element;i++)
1668 for (
int j=0;j<nDOF_trial_element;j++)
1669 elementJacobian_u_u[i][j]=0.0;
1671 for (
int k=0;k<nQuadraturePoints_element;k++)
1675 jac[nSpace*nSpace], jacDet, jacInv[nSpace*nSpace],
1676 u_grad_trial[nDOF_trial_element*nSpace],
1677 dV, u_grad_test_dV[nDOF_test_element*nSpace],
1682 ck.calculateMapping_element(eN,
1686 mesh_trial_ref.data(),
1687 mesh_grad_trial_ref.data(),
1692 ck.calculateH_element(eN,
1694 nodeDiametersArray.data(),
1696 mesh_trial_ref.data(),
1699 dV = fabs(jacDet)*dV_ref[k];
1701 ck.gradTrialFromRef(&u_grad_trial_ref[k*nDOF_trial_element*nSpace],
1705 for (
int j=0;j<nDOF_trial_element;j++)
1706 for (
int I=0;I<nSpace;I++)
1707 u_grad_test_dV[j*nSpace+I] = u_grad_trial[j*nSpace+I]*dV;
1709 epsDiffusion = epsFactDiffusion*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter[eN]);
1710 for(
int i=0;i<nDOF_test_element;i++)
1712 int i_nSpace = i*nSpace;
1713 for(
int j=0;j<nDOF_trial_element;j++)
1715 int j_nSpace = j*nSpace;
1716 elementJacobian_u_u[i][j] +=
1717 ck.NumericalDiffusionJacobian(epsDiffusion,
1718 &u_grad_trial[j_nSpace],
1719 &u_grad_test_dV[i_nSpace]);
1726 for (
int i=0;i<nDOF_test_element;i++)
1728 int eN_i = eN*nDOF_test_element+i;
1729 for (
int j=0;j<nDOF_trial_element;j++)
1731 int eN_i_j = eN_i*nDOF_trial_element+j;
1732 globalJacobian[csrRowIndeces_u_u[eN_i]
1733 + csrColumnOffsets_u_u[eN_i_j]] += elementJacobian_u_u[i][j];
1742 int nQuadraturePoints_elementIn,
1743 int nDOF_mesh_trial_elementIn,
1744 int nDOF_trial_elementIn,
1745 int nDOF_test_elementIn,
1746 int nQuadraturePoints_elementBoundaryIn,
1750 return proteus::chooseAndAllocateDiscretization2D<cppMCorr3P_base,cppMCorr3P,CompKernel>(nSpaceIn,
1751 nQuadraturePoints_elementIn,
1752 nDOF_mesh_trial_elementIn,
1753 nDOF_trial_elementIn,
1754 nDOF_test_elementIn,
1755 nQuadraturePoints_elementBoundaryIn,
1758 return proteus::chooseAndAllocateDiscretization<cppMCorr3P_base,cppMCorr3P,CompKernel>(nSpaceIn,
1759 nQuadraturePoints_elementIn,
1760 nDOF_mesh_trial_elementIn,
1761 nDOF_trial_elementIn,
1762 nDOF_test_elementIn,
1763 nQuadraturePoints_elementBoundaryIn,