Changeset 16231
- Timestamp:
- 09/23/13 10:29:35 (11 years ago)
- Location:
- issm/trunk-jpl/src/c
- Files:
-
- 2 deleted
- 19 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk-jpl/src/c/CMakeLists.txt
r16228 r16231 515 515 ./toolkits/petsc/patches/NewMat.cpp 516 516 ./toolkits/petsc/patches/VecFree.cpp 517 ./toolkits/petsc/patches/PetscMatrixToDoubleMatrix.cpp518 ./toolkits/petsc/patches/PetscVectorToDoubleVector.cpp519 517 ./toolkits/petsc/patches/KSPFree.cpp 520 518 ./toolkits/petsc/patches/MatFree.cpp -
issm/trunk-jpl/src/c/Makefile.am
r16228 r16231 767 767 ./toolkits/petsc/patches/NewMat.cpp\ 768 768 ./toolkits/petsc/patches/VecFree.cpp\ 769 ./toolkits/petsc/patches/PetscMatrixToDoubleMatrix.cpp\770 ./toolkits/petsc/patches/PetscVectorToDoubleVector.cpp\771 769 ./toolkits/petsc/patches/KSPFree.cpp\ 772 770 ./toolkits/petsc/patches/MatFree.cpp\ -
issm/trunk-jpl/src/c/bamg/Curve.cpp
r12821 r16231 20 20 21 21 /*Methods*/ 22 /*FUNCTION Curve::Reverse {{{*/23 void Curve::Reverse() {24 /*reverse the direction of the curve */25 Exchange(FirstEdge,LastEdge);26 Exchange(FirstVertexIndex,LastVertexIndex);27 }28 /*}}}*/29 22 /*FUNCTION Curve::Set {{{*/ 30 23 void Curve::Set(const Curve & rec,const Geometry & Gh ,Geometry & GhNew){ -
issm/trunk-jpl/src/c/bamg/Curve.h
r15067 r16231 19 19 //Methods 20 20 Curve(); 21 void Reverse(void);22 21 void Set(const Curve & rec,const Geometry & Th ,Geometry & ThNew); 23 22 }; -
issm/trunk-jpl/src/c/bamg/GeomEdge.cpp
r15067 r16231 175 175 type |= 64;/*=>6th digit to 1*/ 176 176 }/*}}}*/ 177 /*FUNCTION GeomEdge::Tg{{{*/178 int GeomEdge::Tg(int i) const {179 return i==0 ? TgA() : TgB();180 }/*}}}*/181 177 /*FUNCTION GeomEdge::TgA{{{*/ 182 178 int GeomEdge::TgA() const { -
issm/trunk-jpl/src/c/bamg/GeomEdge.h
r15066 r16231 28 28 R2 F(double theta) const ; // parametrization of the curve edge 29 29 double R1tg(double theta,R2 &t) const ; // 1/radius of curvature + tangente 30 int Tg(int i) const;31 30 int Cracked() const; 32 31 int TgA() const; -
issm/trunk-jpl/src/c/bamg/Mesh.cpp
r16228 r16231 3012 3012 } 3013 3013 /*}}}*/ 3014 /*FUNCTION Mesh::isCracked{{{*/3015 int Mesh::isCracked() const {3016 return NbCrackedVertices != 0;3017 }3018 /*}}}*/3019 3014 /*FUNCTION Mesh::MakeGeomEdgeToEdge{{{*/ 3020 3015 Edge** Mesh::MakeGeomEdgeToEdge() { … … 3817 3812 } 3818 3813 /*}}}*/ 3819 /*FUNCTION Mesh::VerticesRenumber{{{*/3820 void Mesh::VerticesRenumber(long * renu) {3821 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/ReNumberingVertex)*/3822 3823 // warning be carfull because pointer3824 // from on mesh to over mesh3825 // -- so do ReNumbering at the beginning3826 BamgVertex * ve = vertices+nbv;3827 long it,ie,i;3828 3829 _printf_("renumbering triangles\n");3830 for ( it=0;it<nbt;it++)3831 triangles[it].Renumbering(vertices,ve,renu);3832 3833 _printf_("renumbering edges\n");3834 for ( ie=0;ie<nbe;ie++)3835 edges[ie].Renumbering(vertices,ve,renu);3836 3837 _printf_("renumbering vertices on geom\n");3838 for (i=0;i< NbVerticesOnGeomVertex;i++)3839 {3840 BamgVertex *v = VerticesOnGeomVertex[i].meshvertex;3841 if (v>=vertices && v < ve)3842 VerticesOnGeomVertex[i].meshvertex=vertices+renu[GetId(v)];3843 }3844 3845 _printf_("renumbering vertices on edge\n");3846 for (i=0;i< NbVerticesOnGeomEdge;i++)3847 {3848 BamgVertex *v =VerticesOnGeomEdge[i].meshvertex;3849 if (v>=vertices && v < ve)3850 VerticesOnGeomEdge[i].meshvertex=vertices+renu[GetId(v)];3851 }3852 3853 _printf_("renumbering vertices on Bth vertex\n");3854 for (i=0;i< NbVertexOnBThVertex;i++)3855 {3856 BamgVertex *v=VertexOnBThVertex[i].v;3857 if (v>=vertices && v < ve)3858 VertexOnBThVertex[i].v=vertices+renu[GetId(v)];3859 }3860 3861 for (i=0;i< NbVertexOnBThEdge;i++)3862 {3863 BamgVertex *v=VertexOnBThEdge[i].v;3864 if (v>=vertices && v < ve)3865 VertexOnBThEdge[i].v=vertices+renu[GetId(v)];3866 }3867 3868 // move the Vertices without a copy of the array3869 // be carefull not trivial code3870 long j;3871 for ( it=0;it<nbv;it++) // for all sub cycles of the permutation renu3872 if (renu[it] >= 0) // a new sub cycle3873 {3874 i=it;3875 BamgVertex ti=vertices[i],tj;3876 while ( (j=renu[i]) >= 0){3877 // i is old, and j is new3878 renu[i] = -1-renu[i]; // mark3879 tj = vertices[j]; // save new3880 vertices[j]= ti; // new <- old3881 i=j; // next3882 ti = tj;3883 }3884 }3885 if (quadtree){3886 delete quadtree;3887 quadtree = new BamgQuadtree(this);3888 }3889 3890 for ( it=0;it<nbv;it++) renu[i]= -renu[i]-1;3891 }3892 /*}}}*/3893 3814 /*FUNCTION Mesh::SetIntCoor{{{*/ 3894 3815 void Mesh::SetIntCoor(const char * strfrom) { … … 3951 3872 3952 3873 if (number_of_errors) _error_("Fatal error: some triangles have negative areas, see above"); 3953 }3954 /*}}}*/3955 /*FUNCTION Mesh::ShowRegulaty{{{*/3956 void Mesh::ShowRegulaty() const {3957 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr>*/3958 3959 const double sqrt32=sqrt(3.)*0.5;3960 const double aireKh=sqrt32*0.5;3961 D2 Beq(1,0),Heq(0.5,sqrt32);3962 D2xD2 Br(D2xD2(Beq,Heq).t());3963 D2xD2 B1r(Br.inv());3964 double gammamn=1e100,hmin=1e100;3965 double gammamx=0,hmax=0;3966 double beta=1e100;3967 double beta0=0;3968 double alpha2=0;3969 double area=0,Marea=0;3970 // double cf= double(coefIcoor);3971 // double cf2= 6.*cf*cf;3972 int nt=0;3973 for (int it=0;it<nbt;it++)3974 if ( triangles[it].link){3975 Triangle &K=triangles[it];3976 double area3= Area2((R2) K[0],(R2) K[1],(R2) K[2])/6.;3977 area+= area3;3978 D2xD2 B_Kt(K[0],K[1],K[2]);3979 D2xD2 B_K(B_Kt.t());3980 D2xD2 B1K = Br*B_K.inv();3981 D2xD2 BK = B_K*B1r;3982 D2xD2 B1B1 = B1K.t()*B1K;3983 Metric MK(B1B1.x.x,B1B1.x.y,B1B1.y.y);3984 EigenMetric VMK(MK);3985 alpha2 = Max(alpha2,Max(VMK.lambda1/VMK.lambda2,VMK.lambda2/VMK.lambda1));3986 double betaK=0;3987 3988 for (int j=0;j<3;j++)3989 {3990 double he= Norme2(R2(K[j],K[(j+1)%3]));3991 hmin=Min(hmin,he);3992 hmax=Max(hmax,he);3993 BamgVertex & v=K[j];3994 D2xD2 M((Metric)v);3995 betaK += sqrt(M.det());3996 3997 D2xD2 BMB = BK.t()*M*BK;3998 Metric M1(BMB.x.x,BMB.x.y,BMB.y.y);3999 EigenMetric VM1(M1);4000 gammamn=Min3(gammamn,VM1.lambda1,VM1.lambda2);4001 gammamx=Max3(gammamx,VM1.lambda1,VM1.lambda2);4002 }4003 betaK *= area3;// 1/2 (somme sqrt(det))* area(K)4004 Marea+= betaK;4005 beta=min(beta,betaK);4006 beta0=max(beta0,betaK);4007 }4008 area*=3;4009 gammamn=sqrt(gammamn);4010 gammamx=sqrt(gammamx);4011 _printf_(" Adaptmesh info:\n");4012 _printf_(" number of triangles = " << nt << "\n");4013 _printf_(" hmin = " << hmin << ", hmax=" << hmax << "\n");4014 _printf_(" area = " << area << ", M area = " << Marea << ", M area/( |Khat| nt) = " << Marea/(aireKh*nt) << "\n");4015 _printf_(" infinite-regularity(?): min = " << gammamn << ", max = " << gammamx << "\n");4016 _printf_(" anisomax = " << pow(alpha2,0.5) << ", beta max = " << 1./pow(beta/aireKh,0.5) << ", min = " << 1./pow(beta0/aireKh,0.5) << "\n");4017 }4018 /*}}}*/4019 /*FUNCTION Mesh::ShowHistogram{{{*/4020 void Mesh::ShowHistogram() const {4021 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/ShowHistogram)*/4022 4023 const long kmax=10;4024 const double llmin = 0.5,llmax=2;4025 const double lmin=log(llmin),lmax=log(llmax),delta= kmax/(lmax-lmin);4026 long histo[kmax+1];4027 long i,it,k, nbedges =0;4028 for (i=0;i<=kmax;i++) histo[i]=0;4029 for (it=0;it<nbt;it++)4030 if ( triangles[it].link)4031 {4032 4033 for (int j=0;j<3;j++)4034 {4035 Triangle *ta = triangles[it].TriangleAdj(j);4036 if ( !ta || !ta->link || GetId(ta) >= it)4037 {4038 BamgVertex & vP = triangles[it][VerticesOfTriangularEdge[j][0]];4039 BamgVertex & vQ = triangles[it][VerticesOfTriangularEdge[j][1]];4040 if ( !& vP || !&vQ) continue;4041 R2 PQ = vQ.r - vP.r;4042 double l = log(LengthInterpole(vP,vQ,PQ));4043 nbedges++;4044 k = (int) ((l - lmin)*delta);4045 k = Min(Max(k,0L),kmax);4046 histo[k]++;4047 }4048 }4049 }4050 _printf_(" --- Histogram of the unit mesh, nb of edges = " << nbedges << "\n");4051 _printf_(" length of edge in | %% of edge | Nb of edges \n");4052 _printf_(" --------------------+-------------+-------------\n");4053 for (i=0;i<=kmax;i++){4054 if (i==0) _printf_( " " << setw(10) << 0.);4055 else _printf_( " " << setw(10) << exp(lmin+i/delta));4056 if (i==kmax) _printf_(" +inf ");4057 else _printf_( " " << setw(10) << exp(lmin+(i+1)/delta));4058 _printf_("| " << setw(10) << (long((10000. * histo[i])/ nbedges)/100.) << " |\n");4059 _printf_(" " << histo[i] << "\n");4060 }4061 _printf_(" --------------------+-------------+-------------\n");4062 3874 } 4063 3875 /*}}}*/ -
issm/trunk-jpl/src/c/bamg/Mesh.h
r16228 r16231 83 83 long TriangleReferenceList(long*) const; 84 84 void TriangleIntNumbering(long* renumbering); 85 void ShowHistogram() const;86 85 void CrackMesh(BamgOpts* bamgopts); 87 void ShowRegulaty() const;88 86 void SmoothMetric(double raisonmax) ; 89 87 void BoundAnisotropy(double anisomax,double hminaniso= 1e-100) ; … … 96 94 long InsertNewPoints(long nbvold,long & NbTSwap) ; 97 95 void TrianglesRenumberBySubDomain(bool justcompress=false); 98 void VerticesRenumber(long * renu);99 96 void SmoothingVertex(int =3,double=0.3); 100 97 Metric MetricAt (const R2 &) const; … … 118 115 void BuildMetric1(BamgOpts* bamgopts); 119 116 void AddGeometryMetric(BamgOpts* bamgopts); 120 int isCracked() const;121 117 void BuildGeometryFromMesh(BamgOpts* bamgopts=NULL); 122 118 void ReconstructExistingMesh(); -
issm/trunk-jpl/src/c/bamg/Triangle.cpp
r15104 r16231 331 331 void Triangle::SetAllFlag(int a,int f){ 332 332 AdjEdgeIndex[a] = (AdjEdgeIndex[a] &3) + (1020 & f); 333 }/*}}}*/334 /*FUNCTION Triangle::SetDet{{{*/335 void Triangle::SetDet() {336 if(vertices[0] && vertices[1] && vertices[2]) det = bamg::det(*vertices[0],*vertices[1],*vertices[2]);337 else det = -1;338 333 }/*}}}*/ 339 334 /*FUNCTION Triangle::SetHidden{{{*/ -
issm/trunk-jpl/src/c/bamg/Triangle.h
r15066 r16231 63 63 void SetMarkUnSwap(int a); 64 64 void SetUnMarkUnSwap(int a); 65 void SetDet();66 65 67 66 //Inline methods -
issm/trunk-jpl/src/c/classes/Elements/Penta.cpp
r16218 r16231 3029 3029 /*clean-up*/ 3030 3030 delete gauss; 3031 }3032 /*}}}*/3033 /*FUNCTION Penta::ReduceMatrixFS {{{*/3034 void Penta::ReduceMatrixFS(IssmDouble* Ke_reduced, IssmDouble* Ke_temp){3035 3036 int i,j;3037 IssmDouble Kii[24][24];3038 IssmDouble Kib[24][3];3039 IssmDouble Kbb[3][3];3040 IssmDouble Kbi[3][24];3041 IssmDouble Kbbinv[3][3];3042 IssmDouble Kright[24][24];3043 3044 /*Create the four matrices used for reduction */3045 for(i=0;i<24;i++){3046 for(j=0;j<24;j++){3047 Kii[i][j]=*(Ke_temp+27*i+j);3048 }3049 for(j=0;j<3;j++){3050 Kib[i][j]=*(Ke_temp+27*i+j+24);3051 }3052 }3053 for(i=0;i<3;i++){3054 for(j=0;j<24;j++){3055 Kbi[i][j]=*(Ke_temp+27*(i+24)+j);3056 }3057 for(j=0;j<3;j++){3058 Kbb[i][j]=*(Ke_temp+27*(i+24)+j+24);3059 }3060 }3061 3062 /*Inverse the matrix corresponding to bubble part Kbb */3063 Matrix3x3Invert(&Kbbinv[0][0], &Kbb[0][0]);3064 3065 /*Multiply matrices to create the reduce matrix Ke_reduced */3066 TripleMultiply(&Kib[0][0],24,3,0,3067 &Kbbinv[0][0],3,3,0,3068 &Kbi[0][0],3,24,0,3069 &Kright[0][0],0);3070 3071 /*Affect value to the reduced matrix */3072 for(i=0;i<24;i++) for(j=0;j<24;j++) *(Ke_reduced+24*i+j)=Kii[i][j]-Kright[i][j];3073 }3074 /*}}}*/3075 /*FUNCTION Penta::ReduceVectorFS {{{*/3076 void Penta::ReduceVectorFS(IssmDouble* Pe_reduced, IssmDouble* Ke_temp, IssmDouble* Pe_temp){3077 3078 int i,j;3079 IssmDouble Pi[24];3080 IssmDouble Pb[3];3081 IssmDouble Kbb[3][3];3082 IssmDouble Kib[24][3];3083 IssmDouble Kbbinv[3][3];3084 IssmDouble Pright[24];3085 3086 /*Create the four matrices used for reduction */3087 for(i=0;i<24;i++) Pi[i]=*(Pe_temp+i);3088 for(i=0;i<3;i++) Pb[i]=*(Pe_temp+i+24);3089 for(j=0;j<3;j++){3090 for(i=0;i<24;i++){3091 Kib[i][j]=*(Ke_temp+3*i+j);3092 }3093 for(i=0;i<3;i++){3094 Kbb[i][j]=*(Ke_temp+3*(i+24)+j);3095 }3096 }3097 3098 /*Inverse the matrix corresponding to bubble part Kbb */3099 Matrix3x3Invert(&Kbbinv[0][0], &Kbb[0][0]);3100 3101 /*Multiply matrices to create the reduce matrix Ke_reduced */3102 TripleMultiply(&Kib[0][0],24,3,0,3103 &Kbbinv[0][0],3,3,0,3104 &Pb[0],3,1,0,&Pright[0],0);3105 3106 /*Affect value to the reduced matrix */3107 for(i=0;i<24;i++) *(Pe_reduced+i)=Pi[i]-Pright[i];3108 3031 } 3109 3032 /*}}}*/ -
issm/trunk-jpl/src/c/classes/Elements/Penta.h
r16208 r16231 224 224 bool NoIceInElement(void); 225 225 IssmDouble MinEdgeLength(IssmDouble xyz_list[6][3]); 226 void ReduceMatrixFS(IssmDouble* Ke_reduced, IssmDouble* Ke_temp);227 void ReduceVectorFS(IssmDouble* Pe_reduced, IssmDouble* Ke_temp, IssmDouble* Pe_temp);228 226 void SetClone(int* minranks); 229 227 Tria* SpawnTria(int location); -
issm/trunk-jpl/src/c/classes/Loads/Riftfront.cpp
r16164 r16231 733 733 } 734 734 /*}}}*/ 735 /*FUNCTION Riftfront::IsMaterialStable {{{*/736 int Riftfront::IsMaterialStable(void){737 738 IssmDouble converged=0;739 740 this->inputs->GetInputValue(&converged,ConvergedEnum);741 742 if(reCast<int,IssmDouble>(converged)){743 /*ok, material non-linearity has converged. If that was already the case, we keep744 * constraining the rift front. If it was not, and this is the first time the material745 * has converged, we start constraining now!: */746 this->material_converged=1;747 }748 749 return this->material_converged;750 }751 /*}}}*/752 735 /*FUNCTION Riftfront::MaxPenetration {{{*/ 753 736 int Riftfront::MaxPenetration(IssmDouble* ppenetration){ -
issm/trunk-jpl/src/c/classes/Loads/Riftfront.h
r16125 r16231 96 96 int MaxPenetration(IssmDouble* ppenetration); 97 97 int PotentialUnstableConstraint(int* punstable); 98 int IsMaterialStable(void);99 98 bool IsInput(int name); 100 99 /*}}}*/ -
issm/trunk-jpl/src/c/modules/ConstraintsStatex/ConstraintsStateLocal.h
r16206 r16231 21 21 void RiftFreezeConstraints(Loads* loads,int analysis_type); 22 22 23 /*rifts, trial and errors: */24 int RiftIsPreStable(Loads* loads);25 void RiftSetPreStable(Loads* loads);26 void RiftPreConstrain(int* pnum_unstable_constraints,Loads* loads);27 void RiftMaxPenetrationInInputs(Loads* loads);28 int RiftPotentialUnstableConstraints(Loads* loads);29 int RiftIsMaterialStable(Loads* loads);30 31 23 #endif /* _CONSTRAINTSSTATEX_H */ -
issm/trunk-jpl/src/c/modules/ConstraintsStatex/RiftConstraintsState.cpp
r16164 r16231 159 159 } 160 160 /*}}}*/ 161 162 /*diverse trials and errors: */163 /*RiftIsMaterialStable(Loads* loads){{{*/164 int RiftIsMaterialStable(Loads* loads){165 166 int i;167 168 Riftfront* riftfront=NULL;169 int found=0;170 int mpi_found=0;171 172 /*go though loads, and if non-linearity of the material has converged, let all penalties know: */173 for (i=0;i<loads->Size();i++){174 175 if(RiftfrontEnum==loads->GetEnum(i)){176 177 riftfront=(Riftfront*)loads->GetObjectByOffset(i);178 179 if (riftfront->IsMaterialStable()){180 found=1;181 /*do not break! all penalties should get informed the non-linearity converged!*/182 }183 }184 }185 186 ISSM_MPI_Reduce (&found,&mpi_found,1,ISSM_MPI_INT,ISSM_MPI_SUM,0,IssmComm::GetComm() );187 ISSM_MPI_Bcast(&mpi_found,1,ISSM_MPI_INT,0,IssmComm::GetComm());188 found=mpi_found;189 190 return found;191 }192 /*}}}*/193 /*RiftIsPreStable(Loads* loads){{{*/194 int RiftIsPreStable(Loads* loads){195 196 int i;197 198 Riftfront* riftfront=NULL;199 int found=0;200 int mpi_found=0;201 202 /*go though loads, and figure out if one of the penpair loads is still not stable: */203 for (i=0;i<loads->Size();i++){204 205 if(RiftfrontEnum==loads->GetEnum(i)){206 207 riftfront=(Riftfront*)loads->GetObjectByOffset(i);208 209 if (riftfront->PreStable()==0){210 found=1;211 break;212 }213 }214 }215 216 ISSM_MPI_Reduce (&found,&mpi_found,1,ISSM_MPI_INT,ISSM_MPI_SUM,0,IssmComm::GetComm() );217 ISSM_MPI_Bcast(&mpi_found,1,ISSM_MPI_INT,0,IssmComm::GetComm());218 found=mpi_found;219 220 if (found){221 /*We found an unstable constraint. : */222 return 0;223 }224 else{225 return 1;226 }227 }228 /*}}}*/229 /*RiftSetPreStable(Loads* loads){{{*/230 void RiftSetPreStable(Loads* loads){231 232 /*go though loads, and set loads to pre stable.:*/233 for(int i=0;i<loads->Size();i++){234 if(RiftfrontEnum==loads->GetEnum(i)){235 Riftfront* riftfront=(Riftfront*)loads->GetObjectByOffset(i);236 riftfront->SetPreStable();237 }238 }239 }240 /*}}}*/241 /*RiftPreConstrain(int* pnum_unstable_constraints,Loads* loads){{{*/242 void RiftPreConstrain(int* pnum_unstable_constraints,Loads* loads){243 244 int i;245 246 /* generic object pointer: */247 Riftfront* riftfront=NULL;248 249 int unstable;250 int sum_num_unstable_constraints;251 int num_unstable_constraints=0;252 253 /*Enforce constraints: */254 for (i=0;i<loads->Size();i++){255 256 if (RiftfrontEnum==loads->GetEnum(i)){257 258 riftfront=(Riftfront*)loads->GetObjectByOffset(i);259 260 riftfront->PreConstrain(&unstable);261 262 num_unstable_constraints+=unstable;263 }264 }265 266 ISSM_MPI_Reduce (&num_unstable_constraints,&sum_num_unstable_constraints,1,ISSM_MPI_INT,ISSM_MPI_SUM,0,IssmComm::GetComm() );267 ISSM_MPI_Bcast(&sum_num_unstable_constraints,1,ISSM_MPI_INT,0,IssmComm::GetComm());268 num_unstable_constraints=sum_num_unstable_constraints;269 270 /*Assign output pointers: */271 *pnum_unstable_constraints=num_unstable_constraints;272 273 }274 /*}}}*/275 /*RiftMaxPenetrationInInputs(Loads* loads){{{*/276 void RiftMaxPenetrationInInputs(Loads* loads){277 278 int i;279 280 /* generic object pointer: */281 Riftfront* riftfront=NULL;282 283 /*rift penetration: */284 IssmDouble max_penetration=0;285 IssmDouble mpi_max_penetration;286 IssmDouble penetration;287 288 /*Ok, we are going to find the node pairs which are not penetrating, even though they289 * are penalised. We will release only the one with has least <0 penetration. : */290 291 max_penetration=0;292 for (i=0;i<loads->Size();i++){293 294 if (RiftfrontEnum==loads->GetEnum(i)){295 296 riftfront=(Riftfront*)loads->GetObjectByOffset(i);297 298 riftfront->MaxPenetration(&penetration);299 300 if (penetration>max_penetration)max_penetration=penetration;301 }302 }303 304 ISSM_MPI_Reduce (&max_penetration,&mpi_max_penetration,1,ISSM_MPI_DOUBLE,ISSM_MPI_MAX,0,IssmComm::GetComm() );305 ISSM_MPI_Bcast(&mpi_max_penetration,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());306 max_penetration=mpi_max_penetration;307 308 /*feed max_penetration to inputs: */309 for(i=0;i<loads->Size();i++){310 Load* load=(Load*)loads->GetObjectByOffset(i);311 load->InputUpdateFromVector(&max_penetration,MaxPenetrationEnum,ConstantEnum);312 }313 }314 /*}}}*/315 /*RiftPotentialUnstableConstraints(Loads* loads){{{*/316 int RiftPotentialUnstableConstraints(Loads* loads){317 318 int i;319 320 /* generic object pointer: */321 Riftfront* riftfront=NULL;322 323 /*Ok, we are going to find the node pairs which are not penetrating, even though they324 * are penalised. We will release only the one with has least <0 penetration. : */325 int unstable=0;326 int sum_num_unstable_constraints=0;327 int num_unstable_constraints=0;328 329 for (i=0;i<loads->Size();i++){330 331 if (RiftfrontEnum==loads->GetEnum(i)){332 333 riftfront=(Riftfront*)loads->GetObjectByOffset(i);334 335 riftfront->PotentialUnstableConstraint(&unstable);336 337 num_unstable_constraints+=unstable;338 }339 }340 341 ISSM_MPI_Reduce (&num_unstable_constraints,&sum_num_unstable_constraints,1,ISSM_MPI_INT,ISSM_MPI_SUM,0,IssmComm::GetComm() );342 ISSM_MPI_Bcast(&sum_num_unstable_constraints,1,ISSM_MPI_INT,0,IssmComm::GetComm());343 num_unstable_constraints=sum_num_unstable_constraints;344 345 return num_unstable_constraints;346 }347 /*}}}*/ -
issm/trunk-jpl/src/c/shared/TriMesh/TriMeshUtils.cpp
r16218 r16231 490 490 return noerr; 491 491 }/*}}}*/ 492 /*FUNCTION RemoveRifts{{{*/493 int RemoveRifts(int** pindex,double** px,double** py,int* pnods,int** psegments,int* pnumsegs,int numrifts1,int* rifts1numsegs,int** rifts1segments,int** rifts1pairs,int nel){494 495 int noerr=1;496 int i,j,k,counter,counter1,counter2;497 498 /*intermediary: */499 int *riftsegments = NULL;500 int *riftpairs = NULL;501 int node1,node2,node3,node4,temp_node;502 int el2;503 int newnods; //temporary # node counter.504 double xmin,ymin;505 double *xreal = NULL;506 double *yreal = NULL;507 int *nodes = NULL;508 int *mergingnodes = NULL;509 int max_size;510 int redundant;511 512 /*Recover input: */513 int *index = *pindex;514 double *x = *px;515 double *y = *py;516 int nods = *pnods; ;517 int *segments = *psegments;518 int numsegs = *pnumsegs;519 520 /*initialize newnods : */521 newnods=nods;522 523 /*Figure out a unique value to flag x and y for node removal: */524 xmin=x[0];525 ymin=y[0];526 for (i=0;i<nods;i++){527 if (x[i]<xmin)xmin=x[i];528 if (y[i]<ymin)ymin=y[i];529 }530 xmin=xmin-fabs(xmin);531 ymin=ymin-fabs(ymin);532 533 /*Initialize two arrays, one for nodes that are going to be merged, the other with corresponding nodes being merge into: */534 max_size=0;535 for (i=0;i<numrifts1;i++){536 max_size+=rifts1numsegs[i];537 }538 nodes=xNew<int>(max_size);539 mergingnodes=xNew<int>(max_size);540 541 /*Go through the rifts segments, and identify which node we are going to merge with its counterpart on the other side542 *of the rift. The way we identify this node is by looking at the element pairs, and the corresponding nodes: */543 counter=0;544 for (i=0;i<numrifts1;i++){545 riftsegments=rifts1segments[i];546 riftpairs=rifts1pairs[i];547 for (j=0;j<rifts1numsegs[i];j++){548 el2=riftpairs[2*j+1];549 node1=(int)*(riftsegments+3*j+0);550 node2=(int)*(riftsegments+3*j+1);551 /*Summary, el1 and el2 are facing one another across the rift. node1 and node2 belong to el1 and552 *are located on the rift. Find node3 and node4, nodes belonging to el2 and located on the rift: */553 for (k=0;k<rifts1numsegs[i];k++){554 if (*(riftsegments+3*k+2)==el2){555 node3=*(riftsegments+3*k+0);556 node4=*(riftsegments+3*k+1);557 break;558 }559 }560 /* Make sure node3 faces node1 and node4 faces node2: */561 if ((x[node1]==x[node4]) && (y[node1]==y[node4])){562 /*Swap node3 and node4:*/563 temp_node=node3;564 node3=node4;565 node4=temp_node;566 }567 /* Is any of these two node pairs on the tip of a rift, in which case, we don't include it in nodes or mergingnodes: */568 if ((node1==node3) || (node2==node4)){569 if(node1!=node3){570 /*Add node1 and node3 to nodes and mergingnodes if they have not already been added: */571 redundant=0;572 for (k=0;k<counter;k++){573 if ((mergingnodes[k]==node1) || (nodes[k]==node1))redundant=1;574 }575 if(!redundant){576 /*Ok, add node1 to nodes, and node3 to mergingnodes: */577 nodes[counter]=node1;578 mergingnodes[counter]=node3;579 counter++;580 }581 }582 if(node2!=node4){583 /*Add node2 and node4 to nodes and mergingnodes if they have not already been added: */584 redundant=0;585 for (k=0;k<counter;k++){586 if ((mergingnodes[k]==node2) || (nodes[k]==node2))redundant=1;587 }588 if(!redundant){589 /*Ok, add node2 to nodes, and node4 to mergingnodes: */590 nodes[counter]=node2;591 mergingnodes[counter]=node4;592 counter++;593 }594 }595 }596 else{597 /*Check that node1 is not already present in the mergingnodes: */598 redundant=0;599 for (k=0;k<counter;k++){600 if ((mergingnodes[k]==node1) || (nodes[k]==node1))redundant=1;601 }602 if(!redundant){603 /*Ok, add node1 to nodes, and node3 to mergingnodes: */604 nodes[counter]=node1;605 mergingnodes[counter]=node3;606 counter++;607 }608 /*Check that node2 is not already present in the mergingnodes: */609 redundant=0;610 for (k=0;k<counter;k++){611 if ((mergingnodes[k]==node1) || (nodes[k]==node1))redundant=1;612 }613 if(!redundant){614 /*Ok, add node2 to nodes, and node4 to mergingnodes: */615 nodes[counter]=node2;616 mergingnodes[counter]=node4;617 counter++;618 }619 }620 }621 }622 623 /*Ok, we have counter pairs of nodes (nodes and mergingnodes): start merging nodes in the triangulation: */624 newnods=nods;625 for (i=0;i<counter;i++){626 node1=nodes[i];627 node3=mergingnodes[i];628 /*Merge node3 into node1: */629 x[node3]=xmin; //flag for later removal from x630 y[node3]=ymin; //flag for later removal from y631 newnods--;632 for (k=0;k<nel;k++){633 if (*(index+3*k+0)==node3)*(index+3*k+0)=node1;634 if (*(index+3*k+1)==node3)*(index+3*k+1)=node1;635 if (*(index+3*k+2)==node3)*(index+3*k+2)=node1;636 }637 }638 639 /*Reallocate x and y: */640 xreal=xReNew<double>(x,nods,newnods);641 yreal=xReNew<double>(y,nods,newnods);642 counter1=0;643 counter2=0;644 for (i=0;i<nods;i++){645 if (x[i]!=xmin){646 xreal[counter1]=x[i];647 counter1++;648 }649 if (y[i]!=ymin){650 yreal[counter2]=y[i];651 counter2++;652 }653 }654 xDelete<double>(x); x=xreal;655 xDelete<double>(y); y=yreal;656 657 /*Assign output pointers:*/658 *pindex=index;659 *px=x;660 *py=y;661 *pnods=newnods;662 *psegments=segments;663 *pnumsegs=numsegs;664 return noerr;665 }/*}}}*/666 492 /*FUNCTION IsRiftPresent{{{*/ 667 493 int IsRiftPresent(int* priftflag,int* pnumrifts,int* segmentmarkerlist,int nsegs){ -
issm/trunk-jpl/src/c/shared/TriMesh/trimesh.h
r14222 r16231 23 23 int UpdateSegments(int** psegments,int** psegmentmarkerlist, int* pnsegs,int* index, double* x,double* y,int* riftsegments,int nriftsegs,int nods,int nel); 24 24 int FindElement(double A,double B,int* index,int nel); 25 int RemoveRifts(int** pindex,double** px,double** py,int* pnods,int** psegments,int* pnumsegs,int numrifts1,int* rifts1numsegs,int** rifts1segments,double** rifts1pairs,int nel);26 25 int IsRiftPresent(int* priftflag,int* pnumrifts,int* segmentmarkerlist,int nsegs); 27 26 int SplitRiftSegments(int** psegments,int** psegmentmarkerlist, int* pnumsegs, int* pnumrifts,int** priftsnumsegs,int*** priftssegments,int numrifts,int nods,int nels); -
issm/trunk-jpl/src/c/toolkits/petsc/patches/petscpatches.h
r16228 r16231 39 39 MatType ISSMToPetscMatrixType(MatrixType type); 40 40 41 void PetscMatrixToDoubleMatrix(double** pmatrix, int* prows, int* pcols,Mat matrix);42 void PetscVectorToDoubleVector(double** pvector, int* prows, Vec vector);43 44 41 #endif
Note:
See TracChangeset
for help on using the changeset viewer.