![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 2917
- Timestamp:
- 01/27/10 10:25:37 (15 years ago)
- Location:
- issm/trunk/src/c/Bamgx
- Files:
-
- 3 edited
-
Metric.h (modified) (1 diff)
-
objects/MatVVP2x2.cpp (modified) (1 diff)
-
objects/Triangles.cpp (modified) (9 diffs)
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk/src/c/Bamgx/Metric.h
r2899 r2917 108 108 109 109 inline void MatVVP2x2::BoundAniso2(const Real8 coef){ 110 if (coef<=1.00000000001) 111 if (lambda1 < lambda2) 112 lambda1 = bamg::Max(lambda1,lambda2*coef); 113 else 114 lambda2 = bamg::Max(lambda2,lambda1*coef); 115 else // a verifier 116 if (lambda1 > lambda2) 117 lambda1 = bamg::Min(lambda1,lambda2*coef); 118 else 119 lambda2 = bamg::Min(lambda2,lambda1*coef); 120 } 110 if (coef<=1.00000000001){ 111 if (lambda1 < lambda2) 112 lambda1 = bamg::Max(lambda1,lambda2*coef); 113 else 114 lambda2 = bamg::Max(lambda2,lambda1*coef); 115 } 116 else{ //TO BE CHECKED 117 if (lambda1 > lambda2) 118 lambda1 = bamg::Min(lambda1,lambda2*coef); 119 else 120 lambda2 = bamg::Min(lambda2,lambda1*coef); 121 } 122 } 121 123 122 124 void ReductionSimultanee( MetricAnIso M1, MetricAnIso M2,double & l1,double & l2, D2xD2 & V) ; -
issm/trunk/src/c/Bamgx/objects/MatVVP2x2.cpp
r2899 r2917 24 24 double c11 = a11*a11, c22 = a22*a22, c21= a21*a21; 25 25 double b=-a11-a22,c=-c21+a11*a22; 26 double 26 double delta = b*b - 4 * c ; 27 27 double n2=(c11+c22+c21); 28 if ( n2 < 1e-30) 29 lambda1=lambda2=0,v.x=1,v.y=0; 30 else if (delta < eps*n2) 31 { 28 29 //if norm(M)<10^30 -> M=zeros(2,2) 30 if ( n2 < 1e-30) lambda1=lambda2=0,v.x=1,v.y=0; 31 32 //if ??? 33 else if (delta < eps*n2){ 32 34 lambda1=lambda2=-b/2, v.x=1,v.y=0; 33 } 34 else 35 { // --- construction de 2 vecteur dans (Im ( A - D(i) Id) ortogonal 36 delta = sqrt(delta); 37 lambda1 = (-b-delta)/2.0,lambda2 = (-b+delta)/2.0; 35 } 36 37 //general case: construction of 2 eigen vectors 38 else { 39 delta = sqrt(delta); 40 lambda1 = (-b-delta)/2.0,lambda2 = (-b+delta)/2.0; 38 41 double v0 = a11-lambda1, v1 = a21,v2 = a22 - lambda1; 39 42 double s0 = v0*v0 + v1*v1, s1 = v1*v1 +v2*v2; 40 41 43 if(s1 < s0) 42 44 s0=sqrt(s0),v.x=v1/s0,v.y=-v0/s0; 43 45 else 44 46 s1=sqrt(s1),v.x=v2/s1,v.y=-v1/s1; 45 47 }; 46 48 } 47 49 /*}}}1*/ -
issm/trunk/src/c/Bamgx/objects/Triangles.cpp
r2914 r2917 854 854 /*FUNCTION Triangles::BoundAnisotropy{{{1*/ 855 855 void Triangles::BoundAnisotropy(Real8 anisomax,Real8 hminaniso) { 856 856 857 long int verbosity=0; 857 858 858 double lminaniso = 1/ (Max(hminaniso*hminaniso,1e-100)); 859 860 //display info 859 861 if (verbosity > 1) printf(" BoundAnisotropy by %g\n",anisomax); 860 Real8 h1=1.e30,h2=1e-30,rx=0; 862 863 Real8 h1=1.e30,h2=1e-30; 861 864 Real8 coef = 1./(anisomax*anisomax); 862 Real8 hn1=1.e30,hn2=1e-30,rnx =1.e-30 ;863 for (Int4 i=0;i<nbv;i++) 864 {865 865 Real8 hn1=1.e30,hn2=1e-30,rnx =1.e-30,rx=0; 866 867 //loop over all vertices 868 for (Int4 i=0;i<nbv;i++){ 866 869 MatVVP2x2 Vp(vertices[i]); 867 870 double lmax=Vp.lmax(); 868 h1=Min(h1,Vp.lmin()); 869 h2=Max(h2,Vp.lmax()); 870 rx = Max(rx,Vp.Aniso2()); 871 872 Vp *= Min(lminaniso,lmax)/lmax; 873 871 h1 =Min(h1,Vp.lmin()); 872 h2 =Max(h2,Vp.lmax()); 873 Vp*=Min(lminaniso,lmax)/lmax; 874 874 Vp.BoundAniso2(coef); 875 876 hn1=Min(hn1,Vp.lmin());877 hn2=Max(hn2,Vp.lmax());878 rnx = Max(rnx,Vp.Aniso2());879 880 881 875 vertices[i].m = Vp; 882 876 883 } 884 877 //info to be displayed 878 if (verbosity>2){ 879 hn1=Min(hn1,Vp.lmin()); 880 hn2=Max(hn2,Vp.lmax()); 881 rx =Max(rx,Vp.Aniso2()); 882 rnx= Max(rnx,Vp.Aniso2()); 883 } 884 } 885 886 //display info 885 887 if (verbosity>2){ 886 888 printf(" input: Hmin = %g, Hmax = %g, factor of anisotropy max = %g\n",pow(h2,-0.5),pow(h1,-0.5),pow(rx,0.5)); … … 890 892 /*}}}1*/ 891 893 /*FUNCTION Triangles::ConsGeometry{{{1*/ 892 void Triangles::ConsGeometry(Real8 cutoffradian,int *equiedges) // construct a geometry if no geo 893 { 894 void Triangles::ConsGeometry(Real8 cutoffradian,int *equiedges){ 894 895 // if equiedges existe taille nbe 895 896 // equiedges[i]/2 == i original … … 2729 2730 /*FUNCTION Triangles::IntersectConsMetric{{{1*/ 2730 2731 void Triangles::IntersectConsMetric(BamgOpts* bamgopts){ 2731 /*const double* s,const Int4 nbsol,const int * typsols,2732 * const Real8 hmin1,const Real8 hmax1,const Real8 coef,2733 * const Real8 anisomax ,const Real8 CutOff,const int NbJacobi,2734 * const int Rescaling,const double power,const int Hessiantype)*/2735 2732 // the array of solution s is store 2736 2733 // sol0,sol1,...,soln on vertex 0 2737 2734 // sol0,sol1,...,soln on vertex 1 2738 2735 // etc. 2739 // choise = 0 => H is computed with green formula2740 // otherwise=> H is computed from P2 on 4T2736 // Hessiantype = 0 => H is computed with green formula 2737 // Hessiantype = 1 => H is computed from P2 on 4T 2741 2738 2742 2739 /*Options*/ 2743 2740 const int dim = 2; 2744 2741 int AbsError; 2745 double* s ;2742 double* s=NULL; 2746 2743 Int4 nbsol; 2747 int * typsols;2744 int* typsols=NULL; 2748 2745 Real8 hmin1; 2749 2746 Real8 hmax1; … … 2759 2756 /*Recover options*/ 2760 2757 verbosity=bamgopts->verbose; 2761 s=bamgopts->field;2762 nbsol=1; //for now, only one field2763 typsols=0; // only one dof per node2764 2758 AbsError=bamgopts->AbsError; 2765 2759 CutOff=bamgopts->cutoff; … … 2776 2770 coef=sqrt(bamgopts->err)*coef; 2777 2771 2772 /*Get and process fields*/ 2773 s=bamgopts->field; 2774 nbsol=1; //for now, only one field 2775 typsols=(int*)xmalloc(1*sizeof(int)); 2776 typsols[0]=0; // only one dof per node 2777 2778 //sizeoftype = {1,2,3,4} 2778 2779 int sizeoftype[] = { 1, dim ,dim * (dim+1) / 2, dim * dim } ; 2779 2780 2780 // computation of the n b of field2781 // computation of the number of fields 2781 2782 Int4 ntmp = 0; 2782 2783 if (typsols){ 2784 //if there is more than one dof per vertex for one field 2785 //increase ntmp to take into account all the fields 2786 //if nbsol=1 2783 2787 for (Int4 i=0;i<nbsol;i++) ntmp += sizeoftype[typsols[i]]; 2784 2788 } 2785 2789 else ntmp = nbsol; 2786 2790 2787 // n is the total number of fields 2788 2791 //n is the total number of fields 2789 2792 const Int4 n = ntmp; 2790 2793 2791 Int4 i,k,iA,iB,iC,iv; 2792 R2 O(0,0); 2793 int RelativeMetric = CutOff>1e-30; 2794 Real8 hmin = Max(hmin1,MinimalHmin()); 2795 Real8 hmax = Min(hmax1,MaximalHmax()); 2796 Real8 coef2 = 1/(coef*coef); 2797 2794 //initialization of some variables 2795 Int4 i,k,iA,iB,iC,iv; 2796 R2 O(0,0); 2797 int RelativeMetric = CutOff>1e-30; 2798 Real8 hmin = Max(hmin1,MinimalHmin()); 2799 Real8 hmax = Min(hmax1,MaximalHmax()); 2800 Real8 coef2 = 1/(coef*coef); 2801 double* ss=(double*)s; 2802 double sA,sB,sC; 2803 Real8* detT = new Real8[nbt]; 2804 Real8* Mmass= new Real8[nbv]; 2805 Real8* Mmassxx= new Real8[nbv]; 2806 Real8* dxdx= new Real8[nbv]; 2807 Real8* dxdy= new Real8[nbv]; 2808 Real8* dydy= new Real8[nbv]; 2809 Real8* workT= new Real8[nbt]; 2810 Real8* workV= new Real8[nbv]; 2811 int* OnBoundary = new int[nbv]; 2812 2813 //display infos 2798 2814 if(verbosity>1) { 2799 2815 printf(" Construction of Metric: number of field: %i (nbt=%i, nbv=%i)\n",n,nbt,nbv); 2800 2816 printf(" coef = %g\n",coef); 2801 2817 printf(" hmin = %g hmax = %g\n",hmin,hmax); 2802 printf(" anisomax = %g nb Jacobi = %i, power = % i\n",anisomax,NbJacobi,power);2818 printf(" anisomax = %g nb Jacobi = %i, power = %g\n",anisomax,NbJacobi,power); 2803 2819 if (RelativeMetric) printf(" RelativeErr with CutOff= %g\n",CutOff); 2804 2820 else printf(" Absolute error\n"); 2805 2821 } 2806 double* ss=(double*)s; 2807 2808 double sA,sB,sC; 2809 2810 Real8 *detT = new Real8[nbt]; 2811 Real8 *Mmass= new Real8[nbv]; 2812 Real8 *Mmassxx= new Real8[nbv]; 2813 Real8 *dxdx= new Real8[nbv]; 2814 Real8 *dxdy= new Real8[nbv]; 2815 Real8 *dydy= new Real8[nbv]; 2816 Real8 *workT= new Real8[nbt]; 2817 Real8 *workV= new Real8[nbv]; 2818 int *OnBoundary = new int[nbv]; 2822 2823 //initialize Mmass, OnBoundary and Massxx by zero 2819 2824 for (iv=0;iv<nbv;iv++){ 2820 2825 Mmass[iv]=0; … … 2823 2828 } 2824 2829 2825 for (i=0;i<nbt;i++) 2826 if(triangles[i].link){ // the real triangles 2827 const Triangle &t=triangles[i]; 2828 // coor of 3 vertices 2829 R2 A=t[0]; 2830 R2 B=t[1]; 2831 R2 C=t[2]; 2832 2833 2834 // number of the 3 vertices 2835 iA = Number(t[0]); 2836 iB = Number(t[1]); 2837 iC = Number(t[2]); 2838 2839 Real8 dett = bamg::Area2(A,B,C); 2840 detT[i]=dett; 2841 dett /= 6; 2842 2843 // construction of on boundary 2844 int nbb =0; 2845 for(int j=0;j<3;j++) 2846 { 2847 Triangle *ta=t.Adj(j); 2848 if ( ! ta || !ta->link) // no adj triangle => edge on boundary 2849 OnBoundary[Number(t[VerticesOfTriangularEdge[j][0]])]=1, 2850 OnBoundary[Number(t[VerticesOfTriangularEdge[j][1]])]=1, 2851 nbb++; 2852 } 2853 2854 workT[i] = nbb; 2855 Mmass[iA] += dett; 2856 Mmass[iB] += dett; 2857 Mmass[iC] += dett; 2858 2859 if((nbb==0)|| !Hessiantype){ 2860 Mmassxx[iA] += dett; 2861 Mmassxx[iB] += dett; 2862 Mmassxx[iC] += dett; 2863 } 2864 } 2865 else workT[i]=-1; 2866 2830 //Build detT Mmas Mmassxx workT and OnBoundary 2831 for (i=0;i<nbt;i++){ 2832 2833 //lopp over the real triangles (no boundary elements) 2834 if(triangles[i].link){ 2835 2836 //get current triangle t 2837 const Triangle &t=triangles[i]; 2838 2839 // coor of 3 vertices 2840 R2 A=t[0]; 2841 R2 B=t[1]; 2842 R2 C=t[2]; 2843 2844 // number of the 3 vertices 2845 iA = Number(t[0]); 2846 iB = Number(t[1]); 2847 iC = Number(t[2]); 2848 2849 //compute triangle determinant (2*Area) 2850 Real8 dett = bamg::Area2(A,B,C); 2851 detT[i]=dett; 2852 dett /= 6; 2853 2854 // construction of OnBoundary (flag=1 if on boundary, else 0) 2855 int nbb=0; 2856 for(int j=0;j<3;j++){ 2857 //get adjacent triangle 2858 Triangle *ta=t.Adj(j); 2859 //if there is no adjacent triangle, the edge of the triangle t is on boundary 2860 if ( !ta || !ta->link){ 2861 //mark the two vertices of the edge as OnBoundary 2862 OnBoundary[Number(t[VerticesOfTriangularEdge[j][0]])]=1; 2863 OnBoundary[Number(t[VerticesOfTriangularEdge[j][1]])]=1; 2864 nbb++; 2865 } 2866 } 2867 2868 //number of vertices on boundary for current triangle t 2869 workT[i] = nbb; 2870 2871 //Build Mmass Mmass[i] = Mmass[i] + Area/3 2872 Mmass[iA] += dett; 2873 Mmass[iB] += dett; 2874 Mmass[iC] += dett; 2875 2876 //Build Massxx = Mmass 2877 if((nbb==0) || (Hessiantype==0)){ 2878 Mmassxx[iA] += dett; 2879 Mmassxx[iB] += dett; 2880 Mmassxx[iC] += dett; 2881 } 2882 } 2883 2884 //else: the triangle is a boundary triangle -> workT=-1 2885 else workT[i]=-1; 2886 } 2887 2888 //for all Solution 2867 2889 for (Int4 nusol=0;nusol<nbsol;nusol++) { 2868 //for all Solution2869 2890 2870 2891 Real8 smin=ss[0],smax=ss[0]; 2871 2872 2892 Real8 h1=1.e30,h2=1e-30,rx=0; 2873 2893 Real8 coef = 1./(anisomax*anisomax); 2874 2894 Real8 hn1=1.e30,hn2=1e-30,rnx =1.e-30; 2875 int nbfield = typsols? sizeoftype[typsols[nusol]] : 1; 2876 if (nbfield == 1) 2877 for ( iv=0,k=0; iv<nbv; iv++,k+=n ){ 2878 dxdx[iv]=dxdy[iv]=dydy[iv]=0; 2879 smin=Min(smin,ss[k]); 2880 smax=Max(smax,ss[k]); 2881 } 2895 int nbfield=typsols?sizeoftype[typsols[nusol]]:1; 2896 2897 //only one field 2898 if (nbfield == 1) { 2899 //get min(s), max(s) and initialize Hessian (dxdx,dxdy,dydy) 2900 for ( iv=0,k=0; iv<nbv; iv++,k+=n ){ 2901 dxdx[iv]=dxdy[iv]=dydy[iv]=0; 2902 smin=Min(smin,ss[k]); 2903 smax=Max(smax,ss[k]); 2904 } 2905 } 2906 2907 //vectorial case 2882 2908 else{ 2883 // cas vectoriel2884 for ( iv=0,k=0; iv<nbv; iv++,k+=n ){2909 for (iv=0,k=0;iv<nbv;iv++,k+=n ){ 2910 //compute v = √sum(s[i]^2) 2885 2911 double v=0; 2886 for (int i=0;i<nbfield;i++) 2887 v += ss[k+i]*ss[k+i]; 2912 for (int i=0;i<nbfield;i++){ 2913 v += ss[k+i]*ss[k+i]; 2914 } 2888 2915 v = sqrt(v); 2889 2916 smin=Min(smin,v); … … 2891 2918 } 2892 2919 } 2893 Real8 sdelta =smax-smin;2920 Real8 sdelta=smax-smin; 2894 2921 Real8 absmax=Max(Abs(smin),Abs(smax)); 2895 Real8 cnorm = Rescaling ? coef2/sdelta : coef2; 2896 2922 Real8 cnorm =Rescaling ? coef2/sdelta : coef2; 2923 2924 //display info 2897 2925 if(verbosity>2) printf(" Solution %i, Min = %g, Max = %g, Delta = %g, cnorm = %g, number of fields = %i\n",nusol,smin,smax,sdelta,cnorm,nbfield); 2898 2926 2899 if ( sdelta < 1.0e-10*Max(absmax,1e-20) && (nbfield ==1)) { 2927 //skip constant field 2928 if (sdelta < 1.0e-10*Max(absmax,1e-20) && (nbfield ==1)){ 2900 2929 if (verbosity>2) printf(" Solution %i is constant, skipping...\n"); 2901 2930 continue; 2902 2931 } 2903 2932 2904 double *sf = ss; 2933 //pointer toward ss that is also a pointer toward s (solutions) 2934 double* sf=ss; 2935 2936 //loop over all the fields of the solution 2905 2937 for (Int4 nufield=0;nufield<nbfield;nufield++,ss++){ 2906 for ( iv=0,k=0; iv<nbv; iv++,k+=n ) 2907 dxdx[iv]=dxdy[iv]=dydy[iv]=0; 2908 for (i=0;i<nbt;i++) 2909 if(triangles[i].link) 2910 {// for real all triangles 2911 // coor of 3 vertices 2912 R2 A=triangles[i][0]; 2913 R2 B=triangles[i][1]; 2914 R2 C=triangles[i][2]; 2915 2916 2917 // warning the normal is internal and the 2918 // size is the length of the edge 2919 R2 nAB = Orthogonal(B-A); 2920 R2 nBC = Orthogonal(C-B); 2921 R2 nCA = Orthogonal(A-C); 2922 // remark : nAB + nBC + nCA == 0 2923 2924 // number of the 3 vertices 2925 iA = Number(triangles[i][0]); 2926 iB = Number(triangles[i][1]); 2927 iC = Number(triangles[i][2]); 2928 2929 // for the test of boundary edge 2930 // the 3 adj triangles 2931 Triangle *tBC = triangles[i].TriangleAdj(OppositeEdge[0]); 2932 Triangle *tCA = triangles[i].TriangleAdj(OppositeEdge[1]); 2933 Triangle *tAB = triangles[i].TriangleAdj(OppositeEdge[2]); 2934 2935 // value of the P1 fonction on 3 vertices 2936 sA = ss[iA*n]; 2937 sB = ss[iB*n]; 2938 sC = ss[iC*n]; 2939 2940 R2 Grads = (nAB * sC + nBC * sA + nCA * sB ) /detT[i] ; 2941 if(Hessiantype){ 2942 int nbb = 0; 2943 Real8 dd = detT[i]; 2944 Real8 lla,llb,llc,llf; 2945 Real8 taa[3][3],bb[3]; 2946 // construction of the trans of lin system 2947 for (int j=0;j<3;j++) 2948 { 2949 int ie = OppositeEdge[j]; 2950 TriangleAdjacent ta = triangles[i].Adj(ie); 2951 Triangle *tt = ta; 2952 if (tt && tt->link) 2953 { 2954 Vertex &v = *ta.OppositeVertex(); 2955 R2 V = v; 2956 Int4 iV = Number(v); 2957 Real8 lA = bamg::Area2(V,B,C)/dd; 2958 Real8 lB = bamg::Area2(A,V,C)/dd; 2959 Real8 lC = bamg::Area2(A,B,V)/dd; 2960 taa[0][j] = lB*lC; 2961 taa[1][j] = lC*lA; 2962 taa[2][j] = lA*lB; 2963 lla = lA,llb=lB,llc=lC,llf=ss[iV*n] ; 2964 2965 bb[j] = ss[iV*n] - ( sA*lA + sB*lB + sC*lC ) ; 2938 //ss++ so that for each iteration ss points toward the right field 2939 2940 //initialize the hessian matrix 2941 for ( iv=0,k=0; iv<nbv; iv++,k+=n ) dxdx[iv]=dxdy[iv]=dydy[iv]=0; 2942 2943 //loop over the triangles 2944 for (i=0;i<nbt;i++){ 2945 2946 //for real all triangles 2947 if(triangles[i].link){ 2948 2949 // coor of 3 vertices 2950 R2 A=triangles[i][0]; 2951 R2 B=triangles[i][1]; 2952 R2 C=triangles[i][2]; 2953 2954 //warning: the normal is internal and the size is the length of the edge 2955 R2 nAB = Orthogonal(B-A); 2956 R2 nBC = Orthogonal(C-B); 2957 R2 nCA = Orthogonal(A-C); 2958 //note that : nAB + nBC + nCA == 0 2959 2960 // number of the 3 vertices 2961 iA = Number(triangles[i][0]); 2962 iB = Number(triangles[i][1]); 2963 iC = Number(triangles[i][2]); 2964 2965 // for the test of boundary edge 2966 // the 3 adj triangles 2967 Triangle *tBC = triangles[i].TriangleAdj(OppositeEdge[0]); 2968 Triangle *tCA = triangles[i].TriangleAdj(OppositeEdge[1]); 2969 Triangle *tAB = triangles[i].TriangleAdj(OppositeEdge[2]); 2970 2971 // value of the P1 fonction on 3 vertices 2972 sA = ss[iA*n]; 2973 sB = ss[iB*n]; 2974 sC = ss[iC*n]; 2975 2976 /*The nodal functions are such that for a vertex A: 2977 N_A(x,y)=alphaA x + beta_A y +gamma_A 2978 N_A(A) = 1, N_A(B) = 0, N_A(C) = 0 2979 solving this system of equation (determinant = 2Area(T) != 0 if A,B and C are not inlined) 2980 leads to: 2981 N_A = (xB yC - xC yB + x(yB-yC) +y(xC-xB))/(2*Area(T)) 2982 and this gives: 2983 alpha_A = (yB-yC)/(2*Area(T)) 2984 beta_A = (xC-xB)/(2*Area(T)) 2985 and therefore: 2986 grad N_A = nA / detT 2987 for an interpolation of a solution s: 2988 grad(s) = s * sum_{i=A,B,C} grad(N_i) */ 2989 2990 R2 Grads=(nAB*sC+nBC*sA+nCA*sB)/detT[i]; 2991 2992 //from P2 on 4T to compute Hessian 2993 if(Hessiantype==1){ 2994 int nbb=0; 2995 Real8 dd = detT[i]; 2996 Real8 lla,llb,llc,llf; 2997 Real8 taa[3][3],bb[3]; 2998 2999 // construction of the transpose of lin system 3000 for (int j=0;j<3;j++){ 3001 int ie = OppositeEdge[j]; 3002 TriangleAdjacent ta = triangles[i].Adj(ie); 3003 Triangle* tt = ta; 3004 3005 //if the adjacent triangle is not a boundary triangle: 3006 if (tt && tt->link){ 3007 Vertex &v = *ta.OppositeVertex(); 3008 R2 V = v; 3009 Int4 iV = Number(v); 3010 Real8 lA = bamg::Area2(V,B,C)/dd; 3011 Real8 lB = bamg::Area2(A,V,C)/dd; 3012 Real8 lC = bamg::Area2(A,B,V)/dd; 3013 taa[0][j] = lB*lC; 3014 taa[1][j] = lC*lA; 3015 taa[2][j] = lA*lB; 3016 lla = lA,llb=lB,llc=lC,llf=ss[iV*n] ; 3017 bb[j] = ss[iV*n]-(sA*lA+sB*lB+sC*lC) ; 2966 3018 } 2967 else 2968 { 2969 nbb++; 2970 taa[0][j]=0; 2971 taa[1][j]=0; 2972 taa[2][j]=0; 2973 taa[j][j]=1; 2974 bb[j]=0; 3019 else{ 3020 nbb++; 3021 taa[0][j]=0; 3022 taa[1][j]=0; 3023 taa[2][j]=0; 3024 taa[j][j]=1; 3025 bb[j]=0; 2975 3026 } 2976 3027 } 2977 3028 2978 // resolution of 3x3 lineaire system transpose 2979 Real8 det33 = det3x3(taa[0],taa[1],taa[2]); 2980 Real8 cBC = det3x3(bb,taa[1],taa[2]); 2981 Real8 cCA = det3x3(taa[0],bb,taa[2]); 2982 Real8 cAB = det3x3(taa[0],taa[1],bb); 2983 2984 if (!det33){ 2985 throw ErrorException(__FUNCT__,exprintf("!det33")); 2986 } 2987 // computation of the gradient in the element 2988 2989 // H( li*lj) = grad li grad lj + grad lj grad lj 2990 // grad li = njk / detT ; with i j k =(A,B,C) 2991 Real8 Hxx = cAB * ( nBC.x*nCA.x) + cBC * ( nCA.x*nAB.x) + cCA * (nAB.x*nBC.x); 2992 Real8 Hyy = cAB * ( nBC.y*nCA.y) + cBC * ( nCA.y*nAB.y) + cCA * (nAB.y*nBC.y); 2993 Real8 Hxy = cAB * ( nBC.y*nCA.x) + cBC * ( nCA.y*nAB.x) + cCA * (nAB.y*nBC.x) 2994 + cAB * ( nBC.x*nCA.y) + cBC * ( nCA.x*nAB.y) + cCA * (nAB.x*nBC.y); 2995 Real8 coef = 1.0/(3*dd*det33); 2996 Real8 coef2 = 2*coef; 2997 Hxx *= coef2; 2998 Hyy *= coef2; 2999 Hxy *= coef2; 3000 if(nbb==0) 3001 { 3002 dxdx[iA] += Hxx; 3003 dydy[iA] += Hyy; 3004 dxdy[iA] += Hxy; 3005 3006 dxdx[iB] += Hxx; 3007 dydy[iB] += Hyy; 3008 dxdy[iB] += Hxy; 3009 3010 dxdx[iC] += Hxx; 3011 dydy[iC] += Hyy; 3012 dxdy[iC] += Hxy; 3029 // resolution of 3x3 linear system transpose 3030 Real8 det33 = det3x3(taa[0],taa[1],taa[2]); 3031 Real8 cBC = det3x3(bb,taa[1],taa[2]); 3032 Real8 cCA = det3x3(taa[0],bb,taa[2]); 3033 Real8 cAB = det3x3(taa[0],taa[1],bb); 3034 3035 if (!det33){ 3036 throw ErrorException(__FUNCT__,exprintf("!det33")); 3013 3037 } 3014 3038 // computation of the Hessian in the element 3039 3040 // H( li*lj) = grad li grad lj + grad lj grad lj 3041 // grad li = njk / detT ; with i j k =(A,B,C) 3042 Real8 Hxx = cAB * ( nBC.x*nCA.x) + cBC * ( nCA.x*nAB.x) + cCA * (nAB.x*nBC.x); 3043 Real8 Hyy = cAB * ( nBC.y*nCA.y) + cBC * ( nCA.y*nAB.y) + cCA * (nAB.y*nBC.y); 3044 Real8 Hxy = cAB * ( nBC.y*nCA.x) + cBC * ( nCA.y*nAB.x) + cCA * (nAB.y*nBC.x) 3045 + cAB * ( nBC.x*nCA.y) + cBC * ( nCA.x*nAB.y) + cCA * (nAB.x*nBC.y); 3046 Real8 coef = 1.0/(3*dd*det33); 3047 Real8 coef2 = 2*coef; 3048 Hxx *= coef2; 3049 Hyy *= coef2; 3050 Hxy *= coef2; 3051 if(nbb==0){ 3052 dxdx[iA] += Hxx; 3053 dydy[iA] += Hyy; 3054 dxdy[iA] += Hxy; 3055 3056 dxdx[iB] += Hxx; 3057 dydy[iB] += Hyy; 3058 dxdy[iB] += Hxy; 3059 3060 dxdx[iC] += Hxx; 3061 dydy[iC] += Hyy; 3062 dxdy[iC] += Hxy; 3063 } 3015 3064 } 3016 else { 3017 3018 // if edge on boundary no contribution => normal = 03019 if ( ! tBC || ! tBC->link ) nBC = O; 3020 if ( ! tCA || ! tCA->link ) nCA = O;3021 if ( ! tAB || ! tAB->link ) nAB =O;3022 3023 // remark we forgot a 1/2 because3024 // $\\int_{edge} w_i = 1/2 $ if $i$ is in edge 3025 // 0 if not3026 // if we don't take the boundary3027 // dxdx[iA] += ( nCA.x + nAB.x ) *Grads.x;3028 3029 3030 3031 3032 3033 // warning optimization (1) the divide by 2 is done on the metrixconstruction3034 3035 3036 3037 3038 3039 3040 3065 3066 //Use Green to compute Hessian Matrix 3067 else if (Hessiantype==0){ 3068 3069 // if edge on boundary no contribution => normal = 0 3070 if ( !tBC || !tBC->link ) nBC=O; 3071 if ( !tCA || !tCA->link ) nCA=O; 3072 if ( !tAB || !tAB->link ) nAB=O; 3073 3074 // remark we forgot a 1/2 because 3075 // int_{edge} w_i = 1/2 if i is in edge 3076 // 0 if not 3077 // if we don't take the boundary 3078 dxdx[iA] += ( nCA.x + nAB.x ) *Grads.x; 3079 dxdx[iB] += ( nAB.x + nBC.x ) *Grads.x; 3080 dxdx[iC] += ( nBC.x + nCA.x ) *Grads.x; 3081 3082 //warning optimization (1) the division by 2 is done on the metric construction 3083 dxdy[iA] += (( nCA.y + nAB.y ) *Grads.x + ( nCA.x + nAB.x ) *Grads.y) ; 3084 dxdy[iB] += (( nAB.y + nBC.y ) *Grads.x + ( nAB.x + nBC.x ) *Grads.y) ; 3085 dxdy[iC] += (( nBC.y + nCA.y ) *Grads.x + ( nBC.x + nCA.x ) *Grads.y) ; 3086 3087 dydy[iA] += ( nCA.y + nAB.y ) *Grads.y; 3088 dydy[iB] += ( nAB.y + nBC.y ) *Grads.y; 3089 dydy[iC] += ( nBC.y + nCA.y ) *Grads.y; 3041 3090 } 3042 3091 3043 3092 } // for real all triangles 3093 } 3094 3044 3095 Int4 kk=0; 3045 for ( iv=0,k=0 ; iv<nbv; iv++,k+=n ) 3046 if(Mmassxx[iv]>0) 3047 { 3048 dxdx[iv] /= 2*Mmassxx[iv]; 3049 // warning optimization (1) on term dxdy[iv]*ci/2 3050 dxdy[iv] /= 4*Mmassxx[iv]; 3051 dydy[iv] /= 2*Mmassxx[iv]; 3052 // Compute the matrix with abs(eigen value) 3053 Metric M(dxdx[iv], dxdy[iv], dydy[iv]); 3054 MatVVP2x2 Vp(M); 3055 Vp.Abs(); 3056 M = Vp; 3057 dxdx[iv] = M.a11; 3058 dxdy[iv] = M.a21; 3059 dydy[iv] = M.a22; 3096 for ( iv=0,k=0 ; iv<nbv; iv++,k+=n ){ 3097 if(Mmassxx[iv]>0){ 3098 dxdx[iv] /= 2*Mmassxx[iv]; 3099 // warning optimization (1) on term dxdy[iv]*ci/2 3100 dxdy[iv] /= 4*Mmassxx[iv]; 3101 dydy[iv] /= 2*Mmassxx[iv]; 3102 // Compute the matrix with abs(eigen value) 3103 Metric M(dxdx[iv], dxdy[iv], dydy[iv]); 3104 MatVVP2x2 Vp(M); 3105 Vp.Abs(); 3106 M = Vp; 3107 dxdx[iv] = M.a11; 3108 dxdy[iv] = M.a21; 3109 dydy[iv] = M.a22; 3060 3110 } 3061 else kk++; 3111 else kk++; 3112 } 3062 3113 3063 3114 // correction of second derivative 3064 3115 // by a laplacien 3065 Real8 *d2[3] = {dxdx, dxdy, dydy};3066 Real8 *dd;3116 Real8* d2[3] = {dxdx, dxdy, dydy}; 3117 Real8* dd; 3067 3118 for (int xy = 0;xy<3;xy++) { 3068 3119 dd = d2[xy]; … … 3070 3121 for (int ijacobi=0;ijacobi<Max(NbJacobi,2);ijacobi++){ 3071 3122 for (i=0;i<nbt;i++) 3072 if(triangles[i].link) // the real triangles 3073 { 3123 if(triangles[i].link){// the real triangles 3074 3124 // number of the 3 vertices 3075 3125 iA = Number(triangles[i][0]); … … 3080 3130 cc = Max((Real8) ((Mmassxx[iA]>0)+(Mmassxx[iB]>0)+(Mmassxx[iC]>0)),1.); 3081 3131 workT[i] = (dd[iA]+dd[iB]+dd[iC])/cc; 3132 } 3133 for (iv=0;iv<nbv;iv++) workV[iv]=0; 3134 3135 for (i=0;i<nbt;i++){ 3136 if(triangles[i].link){ // the real triangles 3137 // number of the 3 vertices 3138 iA = Number(triangles[i][0]); 3139 iB = Number(triangles[i][1]); 3140 iC = Number(triangles[i][2]); 3141 Real8 cc = workT[i]*detT[i]; 3142 workV[iA] += cc; 3143 workV[iB] += cc; 3144 workV[iC] += cc; 3082 3145 } 3083 for (iv=0;iv<nbv;iv++) 3084 workV[iv]=0; 3085 3086 for (i=0;i<nbt;i++) 3087 if(triangles[i].link) // the real triangles 3088 { 3089 // number of the 3 vertices 3090 iA = Number(triangles[i][0]); 3091 iB = Number(triangles[i][1]); 3092 iC = Number(triangles[i][2]); 3093 Real8 cc = workT[i]*detT[i]; 3094 workV[iA] += cc; 3095 workV[iB] += cc; 3096 workV[iC] += cc; 3146 } 3147 3148 for (iv=0;iv<nbv;iv++){ 3149 if( ijacobi<NbJacobi || OnBoundary[iv]){ 3150 dd[iv] = workV[iv]/(Mmass[iv]*6); 3097 3151 } 3098 3099 for (iv=0;iv<nbv;iv++) 3100 if( ijacobi<NbJacobi || OnBoundary[iv]) 3101 dd[iv] = workV[iv]/(Mmass[iv]*6); 3152 } 3102 3153 } 3103 3154 } 3104 3155 3105 // constuction of the metrixfrom the Hessian dxdx. dxdy,dydy3156 //constuction of the metric from the Hessian dxdx. dxdy,dydy 3106 3157 Real8 rCutOff=CutOff*absmax;// relative cut off 3107 3158 3108 for ( iv=0,k=0 ; iv<nbv; iv++,k+=n ) 3109 { // for all vertices 3159 //loop over the nodes 3160 for ( iv=0,k=0 ; iv<nbv; iv++,k+=n ){ 3161 3110 3162 MetricIso Miso; 3111 3163 Real8 ci ; 3112 if (RelativeMetric) 3113 { // compute the norm of the solution 3164 3165 // compute norm of the solution 3166 if (RelativeMetric){ 3114 3167 double xx =0,*sfk=sf+k; 3115 for (int ifield=0;ifield<nbfield;ifield++,sfk++) 3116 xx += *sfk* *sfk; 3168 for (int ifield=0;ifield<nbfield;ifield++,sfk++){ 3169 xx += *sfk* *sfk; 3170 } 3117 3171 xx=sqrt(xx); 3118 ci = coef2/Max(xx,rCutOff); 3119 } 3120 else ci = cnorm; 3121 3122 Metric Miv(dxdx[iv]*ci, dxdy[iv]*ci, dydy[iv]*ci); 3172 ci=coef2/Max(xx,rCutOff); 3173 } 3174 else ci=cnorm; 3175 3176 //initialize metric Miv with ci*H 3177 Metric Miv(dxdx[iv]*ci, dxdy[iv]*ci, dydy[iv]*ci); 3178 3179 //Get eigen values and vectors of Miv 3123 3180 MatVVP2x2 Vp(Miv); 3124 3181 3182 //move eigen valuse to their absolute values 3125 3183 Vp.Abs(); 3184 3185 //Allpy a power if requested by user 3126 3186 if(power!=1.0) Vp.pow(power); 3127 3187 3188 //get minimum and maximum eigen values 3128 3189 h1=Min(h1,Vp.lmin()); 3129 3190 h2=Max(h2,Vp.lmax()); 3130 3191 3192 //modify eigen values according to hmin and hmax 3131 3193 Vp.Maxh(hmin); 3132 3194 Vp.Minh(hmax); 3133 3195 3196 //multiply eigen values by coef 3197 Vp.BoundAniso2(coef); 3198 3199 //rebuild Metric from Vp 3200 Metric MVp(Vp); 3201 3202 //Apply Metric to vertex 3203 vertices[iv].m.IntersectWith(MVp); 3204 3205 //info to be displayed 3206 //rx = max(lmax/lmin) (anisotropy ratio) 3134 3207 rx = Max(rx,Vp.Aniso2()); 3135 3136 Vp.BoundAniso2(coef);3137 3138 3208 hn1=Min(hn1,Vp.lmin()); 3139 3209 hn2=Max(hn2,Vp.lmax()); 3140 3210 rnx = Max(rnx,Vp.Aniso2()); 3141 3142 Metric MVp(Vp); 3143 vertices[iv].m.IntersectWith(MVp); 3144 //vertices[iv].m=Vp; TEST: work like that 3145 }// for all vertices 3146 //vertices[nbv-1].m.Echo(); 3147 if (verbosity>-1) { 3211 } 3212 3213 //display info 3214 if (verbosity>2) { 3215 3216 3148 3217 printf(" Field %i of solution %i\n",nufield,nusol); 3149 3218 printf(" before bounding : Hmin = %g, Hmax = %g, factor of anisotropy max = %g\n",pow(h2,-0.5), pow(h1,-0.5), pow(rx,0.5)); 3150 3219 printf(" after bounding : Hmin = %g, Hmax = %g, factor of anisotropy max = %g\n",pow(hn2,-0.5),pow(hn1,-0.5),pow(rnx,0.5)); 3151 3220 } 3152 3221 } // end of for all field 3153 3222 }// end for all solution 3154 3223
Note:
See TracChangeset
for help on using the changeset viewer.
