∇-Nabla: Numerical Analysis BAsed LAnguage
∇ Monotone Non-Linear DDFV
These ∇ jobs are taken from the MNLDDFV example: it is a monotone non-linear finite volume method for approximating diffusion operators on general meshes
// ***************************************************************************** // Si on a des quads et que l'on souhaite un Randomly distorted quadrilateral mesh // ***************************************************************************** ∀ nodes void randomDistortedQuads(void) @ -20.0 if (option_rdq and option_quads) { const ℝ α=option_rdq_α; const ℝ ρ1=drand48()+drand48()-1.0; const ℝ ρ2=drand48()+drand48()-1.0; if (coord.x == 0.0 || coord.x == 1.0) continue; if (coord.y == 0.0 || coord.y == 1.0) continue; coord.x+=α*ρ1*Δl; coord.y+=α*ρ2*Δl; } // ***************************************************************************** // * // ***************************************************************************** ∀ nodes void randomDistortedTriangles(void) @ -20.0 if (option_rdq and option_triangles) { const ℝ α=option_rdq_α; const ℝ ρ1=drand48()+drand48()-1.0; const ℝ ρ2=drand48()+drand48()-1.0; if (coord.x == 0.0 || coord.x == 1.0) continue; if (coord.y == 0.0 || coord.y == 1.0) continue; coord.x+=α*ρ1*Δl; coord.y+=α*ρ2*Δl; } // ***************************************************************************** // * // ***************************************************************************** ∀ nodes void stronglyNonConvexQuads(void) @ -20.0 if (option_sncq and option_quads){ const Integer nid=uid+1; const Integer nbNodesPerLine=sqrtl(globalNbNodes); const ℝ θ=M_PI*option_sncq_θ; const ℝ Δborders=Δl*(1.0+cos(θ))/2.0; ℝ Δ=1.5*Δl; // On saute un noeud sur deux if (!(nid%2)) continue; // On saute la dernière colonne if ((nid%nbNodesPerLine)==0) continue; // On saute la première colonne if (((nid-1)%nbNodesPerLine)==0) continue; // On saute la première et dernière ligne if ((((nid-1)/nbNodesPerLine)==0)) continue; if ((((nid-1)/nbNodesPerLine)==(nbNodesPerLine-1))) continue; // A l'avant dernière colonne, on change les valeures if (((nid+1)%nbNodesPerLine)==0) Δ=Δborders; // A l'avant dernière ligne, on change les valeures if ((((nid-1)/nbNodesPerLine)==(nbNodesPerLine-2))) Δ=Δborders; coord.x+=Δ; coord.y+=Δ; } /////////////////////////////////////////////// // ANIsotropic Diffusion Square withOUT Hole // // But with option_θ that must be tied to 0, // // for A and B to be positives. // /////////////////////////////////////////////// ℝ exact_solution_anisotropic_without_hole_and_null_θ(ℝ³ p){ assert(option_𝜕Ω_temperature==0.0); return option_𝜕Ω_temperature+sin(M_PI*p.x)*sin(M_PI*p.y); } ℝ f_anisotropic_without_hole_and_null_θ(ℝ³ p){ assert(option_θ==0.0); return M_PI*M_PI*((option_k+1.0)*sin(M_PI*p.x)*sin(M_PI*p.y) -(option_k-1.0)*cos(M_PI*p.x)*cos(M_PI*p.y)*sin(2.0*option_θ));} /////////////////////////////////////////////// // ANIsotropic Diffusion Square withOUT Hole // // Found a function so that p and f are >=0 // /////////////////////////////////////////////// ℝ exact_solution_anisotropic_without_hole(ℝ³ p){ const ℝ γ=1.0/atan(½); return γ*atan(½-((p.x-½)²+(p.y-½)²)); } ℝ f_anisotropic_without_hole(ℝ³ p){ const ℝ x=p.x; const ℝ y=p.y; const ℝ γ=1.0/atan(½); const ℝ θ=option_θ; const ℝ k=option_k; const ℝ dnm = (1.0+(x²-x+y²-y)²)²; const ℝ num = -16.0*γ*(k-1.0)*(x-½)*(y-½)*(-x+x²+(y-1.0)*y)*cos(θ)*sin(θ) -8.0*γ*(y-½)²*(x²-x+y*(y-1.0))*(k*cos(θ)²+sin(θ)²) -8.0*γ*(x-½)²*(x²-x+y*(y-1.0))*(cos(θ)²+k*sin(θ)²) +2.0*γ*(1.0+(x²-x+(y-1.0)*y)²)*(k*cos(θ)²+sin(θ)²) +2.0*γ*(1.0+(x²-x+(y-1.0)*y)²)*(cos(θ)²+k*sin(θ)²); return num/dnm; } //////////////////////////////////////////// // ANIsotropic Diffusion Square WITH Hole // //////////////////////////////////////////// ℝ g_with_hole(ℝ³ p){ const ℝ θ_hole = option_𝜕Ω_temperature+2.0; const ℝ θ_bord = option_𝜕Ω_temperature; if (p.x==0.0 || p.x== 1.0) return θ_bord; if (p.y==0.0 || p.y== 1.0) return θ_bord; return θ_hole; } // **************************************************************************** // * Fonction 'qui tourne' pour trouver le bon secteur dans le cas convex // **************************************************************************** ℾ convexLoopOnThisNodeToFindPositiveDualFace(ℾ dbg, Node nd, ℝ³x3 kappa, ℝ³ νs, ℝ³ d, ℝ³ *j, ℝ³ *k, ℝ *pj, ℝ *pk, ℝ *α, ℝ *β, int *face_uid, ℾ *face_swap){ foreach nd face{ if (fnd->nbCell()==2){ const Node n0=fnd->node(0); const Node n1=fnd->node(1); const Cell bC = fnd->backCell(); const Cell fC = fnd->frontCell(); const ℝ³ s=½*(coord[n0]+coord[n1]); const ℝ³ p=cell_mass_center[bC]; const ℝ³ q=cell_mass_center[fC]; const ℾ swap=Sin(q-d,s-d)<0.0; *face_swap=swap; *j=swap?p:q; *k=swap?q:p; *pj=swap?cell_θ[bC]:cell_θ[fC]; *pk=swap?cell_θ[fC]:cell_θ[bC]; }else{ const Cell bC = fnd->cell(0); const ℝ³ p=cell_mass_center[bC]; const ℝ³ q=½*(coord[fnd->node(1)]+coord[fnd->node(0)]); const ℝ³ s=½*(p+q); const ℾ swap=Sin(q-d,s-d)<0.0; *face_swap=swap; *j=swap?p:q; *k=swap?q:p; *pj=swap?cell_θ[bC]:g(q); *pk=swap?g(q):cell_θ[bC]; } *α=n(d,*k)⋅(kappa⨂νs); *β=n(*j,d)⋅(kappa⨂νs); if (!(*α>=0.0 && *β>=0.0)) continue; *face_uid=fnd->uniqueId().asInteger(); return true; } return false; } // Pour les faces internes, tout est bien déterminé // Pas de swap à prévoir pour le coté dual non plus ∀ own inner faces void dDual(void) @ 1.0 if (!option_indirect){ const Integer nid0 = 1+node(0)->uniqueId().asInteger(); const Integer nid1 = 1+node(1)->uniqueId().asInteger(); const ℾ dbg=option_debug_dual; ℝ³ j,k,l,m; ℝ pj,pk,pl,pm; ℝ ad,bd,ae,be; int tail_face_uid, head_face_uid; ℾ tail_face_swap, head_face_swap; const ℝ³ d=coord[0]; const ℝ³ e=coord[1]; const ℝ³ p=cell_mass_center[backCell]; const ℝ³ q=cell_mass_center[frontCell]; const ℝ³ νs=n(q,p); const ℾ convexD=geomComputeTriangleAlgebraicArea(d,q,p)>0.0; const ℾ convexE=geomComputeTriangleAlgebraicArea(e,p,q)>0.0; const ℾ convex=convexD&&convexE; if (!convex) fatal("Direct and !convex here"); if (dbg) info()<<"\33[32m[interiorDualFluxesApproximation] Face #"<<uid<<": "<<nid0<<"-"<<nid1<<"\33[m"; { const ℾ okTail = convexLoopOnThisNodeToFindPositiveDualFace(dbg,node(0),κ,νs,d, &j,&k,&pj,&pk,&ad,&bd, &tail_face_uid, &tail_face_swap); const ℾ okHead = convexLoopOnThisNodeToFindPositiveDualFace(dbg,node(1),κ,-νs,e, &l,&m,&pl,&pm,&ae,&be, &head_face_uid, &head_face_swap); const ℝ Ad=geomComputeTriangleArea(d,j,k); //geomComputeTriangleArea(d,j,e)+geomComputeTriangleArea(e,k,d); const ℝ Ae=geomComputeTriangleArea(e,l,m); //geomComputeTriangleArea(e,l,d)+geomComputeTriangleArea(d,m,e); const ℝ μsd_num=ae*pl+be*pm; const ℝ μse_num=ad*pj+bd*pk; const ℝ μs_denum=Ae*(ad*pj+bd*pk)+Ad*(ae*pl+be*pm); const ℾ null=(μs_denum==0.0); const ℝ μsd=null?½:μsd_num; const ℝ μse=null?½:μse_num; const ℝ μsd_denum=null?Ad:μs_denum; const ℝ μse_denum=null?Ae:μs_denum; if(okHead && okTail){ //info()<<"\t\33[33m ok\33[m"; interior_dual_c_sd = ½*(ad+bd)*μsd/μsd_denum; interior_dual_c_se = ½*(ae+be)*μse/μse_denum; interior_dual_c_sl = interior_dual_c_sm = -∞; interior_dual_c_sj = interior_dual_c_sk = -∞; continue; } if (okHead){ // Ici, c'est node[0] qui est on_𝜕Ω, μ=0.0 //info()<<"\t\33[33m !okTail \33[m"; interior_dual_c_lm = true; interior_dual_face_uid=head_face_uid; interior_dual_face_swap=head_face_swap; interior_dual_c_se = ½*(ae+be)/Ae; interior_dual_c_sd = -∞; interior_dual_c_sl = ½*ae/Ae; interior_dual_c_sm = ½*be/Ae; continue; } if (okTail){ // Ici, c'est node[1] qui est on_𝜕Ω, μ=1.0 //info()<<"\t\33[33m !okHead \33[m"; interior_dual_c_jk = true; interior_dual_face_uid=tail_face_uid; interior_dual_face_swap=tail_face_swap; interior_dual_c_sd = ½*(ad+bd)/Ad; interior_dual_c_se = -∞; interior_dual_c_sj = ½*ad/Ad; interior_dual_c_sk = ½*bd/Ad; continue; } fatal("Should not be there!"); } }