Changeset 16231


Ignore:
Timestamp:
09/23/13 10:29:35 (11 years ago)
Author:
Mathieu Morlighem
Message:

DEL: removed legacy code

Location:
issm/trunk-jpl/src/c
Files:
2 deleted
19 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/CMakeLists.txt

    r16228 r16231  
    515515                                        ./toolkits/petsc/patches/NewMat.cpp
    516516                                        ./toolkits/petsc/patches/VecFree.cpp
    517                                         ./toolkits/petsc/patches/PetscMatrixToDoubleMatrix.cpp
    518                                         ./toolkits/petsc/patches/PetscVectorToDoubleVector.cpp
    519517                                        ./toolkits/petsc/patches/KSPFree.cpp
    520518                                        ./toolkits/petsc/patches/MatFree.cpp

  • issm/trunk-jpl/src/c/Makefile.am

    r16228 r16231  
    767767                                        ./toolkits/petsc/patches/NewMat.cpp\
    768768                                        ./toolkits/petsc/patches/VecFree.cpp\
    769                                         ./toolkits/petsc/patches/PetscMatrixToDoubleMatrix.cpp\
    770                                         ./toolkits/petsc/patches/PetscVectorToDoubleVector.cpp\
    771769                                        ./toolkits/petsc/patches/KSPFree.cpp\
    772770                                        ./toolkits/petsc/patches/MatFree.cpp\

  • issm/trunk-jpl/src/c/bamg/Curve.cpp

    r12821 r16231  
    2020
    2121        /*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         /*}}}*/
    2922        /*FUNCTION Curve::Set {{{*/
    3023        void Curve::Set(const Curve & rec,const Geometry & Gh ,Geometry & GhNew){

  • issm/trunk-jpl/src/c/bamg/Curve.h

    r15067 r16231  
    1919                        //Methods
    2020                        Curve();
    21                         void Reverse(void);
    2221                        void Set(const Curve & rec,const Geometry & Th ,Geometry & ThNew);
    2322        };

  • issm/trunk-jpl/src/c/bamg/GeomEdge.cpp

    r15067 r16231  
    175175                type |= 64;/*=>6th digit to 1*/
    176176        }/*}}}*/
    177           /*FUNCTION GeomEdge::Tg{{{*/
    178         int    GeomEdge::Tg(int i) const  {
    179                 return i==0 ? TgA() : TgB();
    180         }/*}}}*/
    181177        /*FUNCTION GeomEdge::TgA{{{*/
    182178        int    GeomEdge::TgA()     const  {

  • issm/trunk-jpl/src/c/bamg/GeomEdge.h

    r15066 r16231  
    2828                        R2     F(double theta) const ; // parametrization of the curve edge
    2929                        double R1tg(double theta,R2 &t) const ; // 1/radius of curvature + tangente
    30                         int    Tg(int i) const;
    3130                        int    Cracked() const;
    3231                        int    TgA()     const;

  • issm/trunk-jpl/src/c/bamg/Mesh.cpp

    r16228 r16231  
    30123012        }
    30133013        /*}}}*/
    3014         /*FUNCTION Mesh::isCracked{{{*/
    3015         int Mesh::isCracked() const {
    3016                 return NbCrackedVertices != 0;
    3017         }
    3018         /*}}}*/
    30193014        /*FUNCTION Mesh::MakeGeomEdgeToEdge{{{*/
    30203015        Edge** Mesh::MakeGeomEdgeToEdge() {
    … …  
    38173812        }
    38183813        /*}}}*/
    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 pointer
    3824                 // from on mesh to over mesh
    3825                 //  --  so do ReNumbering at the beginning
    3826                 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 array
    3869                 // be carefull not trivial code
    3870                 long j;
    3871                 for ( it=0;it<nbv;it++) // for all sub cycles of the permutation renu
    3872                  if (renu[it] >= 0) // a new sub cycle
    3873                         {
    3874                          i=it;
    3875                          BamgVertex ti=vertices[i],tj;
    3876                          while ( (j=renu[i]) >= 0){
    3877                                  // i is old, and j is new
    3878                                  renu[i] = -1-renu[i]; // mark
    3879                                  tj = vertices[j];     // save new
    3880                                  vertices[j]= ti;      // new <- old
    3881                                  i=j;     // next
    3882                                  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         /*}}}*/
    38933814/*FUNCTION Mesh::SetIntCoor{{{*/
    38943815void Mesh::SetIntCoor(const char * strfrom) {
    … …  
    39513872
    39523873        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");
    40623874}
    40633875/*}}}*/

  • issm/trunk-jpl/src/c/bamg/Mesh.h

    r16228 r16231  
    8383                        long TriangleReferenceList(long*) const;
    8484                        void TriangleIntNumbering(long* renumbering);
    85                         void ShowHistogram() const;
    8685                        void CrackMesh(BamgOpts* bamgopts);
    87                         void ShowRegulaty() const;
    8886                        void SmoothMetric(double raisonmax) ;
    8987                        void BoundAnisotropy(double anisomax,double hminaniso= 1e-100) ;
    … …  
    9694                        long InsertNewPoints(long nbvold,long & NbTSwap) ;
    9795                        void TrianglesRenumberBySubDomain(bool justcompress=false);
    98                         void VerticesRenumber(long * renu);
    9996                        void SmoothingVertex(int =3,double=0.3);
    10097                        Metric MetricAt (const R2 &) const;
    … …  
    118115                        void BuildMetric1(BamgOpts* bamgopts);
    119116                        void AddGeometryMetric(BamgOpts* bamgopts);
    120                         int  isCracked() const;
    121117                        void BuildGeometryFromMesh(BamgOpts* bamgopts=NULL);
    122118                        void ReconstructExistingMesh();

  • issm/trunk-jpl/src/c/bamg/Triangle.cpp

    r15104 r16231  
    331331        void   Triangle::SetAllFlag(int a,int f){
    332332                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;
    338333        }/*}}}*/
    339334        /*FUNCTION Triangle::SetHidden{{{*/

  • issm/trunk-jpl/src/c/bamg/Triangle.h

    r15066 r16231  
    6363                        void              SetMarkUnSwap(int a);
    6464                        void              SetUnMarkUnSwap(int a);
    65                         void              SetDet();
    6665
    6766                        //Inline methods

  • issm/trunk-jpl/src/c/classes/Elements/Penta.cpp

    r16218 r16231  
    30293029        /*clean-up*/
    30303030        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];
    31083031}
    31093032/*}}}*/

  • issm/trunk-jpl/src/c/classes/Elements/Penta.h

    r16208 r16231  
    224224                bool           NoIceInElement(void);
    225225                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);
    228226                void             SetClone(int* minranks);
    229227                Tria*            SpawnTria(int location);

  • issm/trunk-jpl/src/c/classes/Loads/Riftfront.cpp

    r16164 r16231  
    733733}
    734734/*}}}*/
    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 keep
    744                  * constraining the rift front. If it was not, and this is the first time the material
    745                  * has converged, we start constraining now!: */
    746                 this->material_converged=1;
    747         }
    748 
    749         return this->material_converged;
    750 }
    751 /*}}}*/
    752735/*FUNCTION Riftfront::MaxPenetration {{{*/
    753736int   Riftfront::MaxPenetration(IssmDouble* ppenetration){

  • issm/trunk-jpl/src/c/classes/Loads/Riftfront.h

    r16125 r16231  
    9696                int   MaxPenetration(IssmDouble* ppenetration);
    9797                int   PotentialUnstableConstraint(int* punstable);
    98                 int   IsMaterialStable(void);
    9998                bool  IsInput(int name);
    10099                /*}}}*/

  • issm/trunk-jpl/src/c/modules/ConstraintsStatex/ConstraintsStateLocal.h

    r16206 r16231  
    2121void RiftFreezeConstraints(Loads* loads,int analysis_type);
    2222
    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 
    3123#endif  /* _CONSTRAINTSSTATEX_H */

  • issm/trunk-jpl/src/c/modules/ConstraintsStatex/RiftConstraintsState.cpp

    r16164 r16231  
    159159}
    160160/*}}}*/
    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 they
    289          * 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 they
    324          * 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  
    490490        return noerr;
    491491}/*}}}*/
    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 side
    542          *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 and
    552                          *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 x
    630                 y[node3]=ymin; //flag  for later removal from y
    631                 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 }/*}}}*/
    666492/*FUNCTION IsRiftPresent{{{*/
    667493int IsRiftPresent(int* priftflag,int* pnumrifts,int* segmentmarkerlist,int nsegs){

  • issm/trunk-jpl/src/c/shared/TriMesh/trimesh.h

    r14222 r16231  
    2323int UpdateSegments(int** psegments,int** psegmentmarkerlist, int* pnsegs,int* index, double* x,double* y,int* riftsegments,int nriftsegs,int nods,int nel);
    2424int 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);
    2625int IsRiftPresent(int* priftflag,int* pnumrifts,int* segmentmarkerlist,int nsegs);
    2726int 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  
    3939MatType ISSMToPetscMatrixType(MatrixType type);
    4040
    41 void PetscMatrixToDoubleMatrix(double** pmatrix, int* prows, int* pcols,Mat matrix);
    42 void PetscVectorToDoubleVector(double** pvector, int* prows, Vec vector);
    43 
    4441#endif
Note: See TracChangeset for help on using the changeset viewer.