source: issm/oecreview/Archive/16133-16554/ISSM-16236-16237.diff

Last change on this file was 16556, checked in by Mathieu Morlighem, 11 years ago

NEW: added Archive/16133-16554
File size: 14.6 KB

  • ../trunk-jpl/src/c/bamg/Geometry.h

     
    5959                        long             GetId(const GeomEdge *t) const;
    6060                        long             GetId(const Curve *c) const;
    6161                        void             UnMarkEdges();
    62                         GeomEdge *ProjectOnCurve(const Edge &,double,BamgVertex &,VertexOnGeom &) const;
    63                         GeomEdge *Containing(const R2 P, GeomEdge *start) const;
     62                        GeomEdge        *ProjectOnCurve(const Edge &,double,BamgVertex &,VertexOnGeom &) const;
    6463                        void             WriteGeometry(BamgGeom *bamggeom, BamgOpts*bamgopts);
    6564        };
    6665

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

     
    251251                        t->AdjEdgeIndex[aat]=a;
    252252                }
    253253        }/*}}}*/
    254         /*FUNCTION Triangle::SetAllFlag{{{*/
    255         void   Triangle::SetAllFlag(int a,int f){
    256                 AdjEdgeIndex[a] = (AdjEdgeIndex[a] &3) + (1020 & f);
    257         }/*}}}*/
    258254        /*FUNCTION Triangle::SetHidden{{{*/
    259255        void Triangle::SetHidden(int a){
    260256                //Get Adjacent Triangle number a

  • ../trunk-jpl/src/c/bamg/Geometry.cpp

     
    462462        long Geometry::GetId(const Curve * c) const  {
    463463                return c - curves;
    464464        }/*}}}*/
    465         /*FUNCTION Geometry::Containing{{{*/
    466         GeomEdge* Geometry::Containing(const R2 P,  GeomEdge * start) const {
    467                 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshGeom.cpp/Contening)*/
    468 
    469                 GeomEdge* on =start,* pon=0;
    470                 // walk with the cos on geometry
    471                 int counter=0;
    472                 while(pon != on){ 
    473                         counter++;
    474                         _assert_(counter<100);
    475                         pon = on;
    476                         R2 A= (*on)[0];
    477                         R2 B= (*on)[1];
    478                         R2 AB = B-A;
    479                         R2 AP = P-A;
    480                         R2 BP = P-B;
    481                         if ( (AB,AP) < 0)
    482                          on = on->Adj[0];
    483                         else if ( (AB,BP)  > 0)
    484                          on = on->Adj[1];
    485                         else
    486                          return on;
    487                 }
    488                 return on;
    489         }
    490         /*}}}*/
    491465        /*FUNCTION Geometry::PostRead{{{*/
    492466        void Geometry::PostRead(){
    493467                /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshGeom.cpp/AfterRead)*/

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

     
    4545                        int               Locked(int a)const;
    4646                        int               Hidden(int a)const;
    4747                        int               GetAllflag(int a);
    48                         void              SetAllFlag(int a,int f);
    4948                        double            QualityQuad(int a,int option=1) const;
    5049                        short             NuEdgeTriangleAdj(int i) const;
    5150                        AdjacentTriangle  Adj(int i) const;

  • ../trunk-jpl/src/c/bamg/Metric.h

     
    7676                        double lmax()   const;
    7777                        double lmin()   const;
    7878                        double Aniso2() const;
    79                         double Aniso()  const;
    80                         void   BoundAniso(const  double c);
    8179                        inline void BoundAniso2(const double coef);
    8280
    8381                        //operators

  • ../trunk-jpl/src/c/bamg/EigenMetric.cpp

     
    9898        void   EigenMetric::Abs(){
    9999                lambda1=bamg::Abs(lambda1),lambda2=bamg::Abs(lambda2);
    100100        }/*}}}*/
    101         /*FUNCTION EigenMetric::Aniso{{{*/
    102         double EigenMetric::Aniso() const  {
    103                 return sqrt( Aniso2());
    104         }/*}}}*/
    105101        /*FUNCTION EigenMetric::Aniso2{{{*/
    106102        double EigenMetric::Aniso2() const  {
    107103                return lmax()/lmin();
    108104        }/*}}}*/
    109         /*FUNCTION EigenMetric::BoundAniso{{{*/
    110         void   EigenMetric::BoundAniso(const double c){
    111                 BoundAniso2(1/(c*c));
    112         }/*}}}*/
    113105        /*FUNCTION EigenMetric::Echo {{{*/
    114106        void EigenMetric::Echo(void){
    115107

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

     
    138138                                triangles[nbt].color = Tho.subdomains[reft[i]].ReferenceNumber;
    139139                                nbt++;           
    140140                          }
    141                   if (kt!=nbt){
    142                           _error_("kt!=nbt");
    143                   }
    144141                  if (nbt==0 && nbv==0) {
    145142                          _error_("All triangles have been removed");
    146143                  }
    147144                  delete [] kk;
    148145                  delete [] reft;
    149146                  delete [] refv;
    150                   //double cutoffradian = 10.0/180.0*Pi;
    151147                  BuildGeometryFromMesh(bamgopts);
    152148                  Gh.PostRead();
    153149                  SetIntCoor();
    154150                  ReconstructExistingMesh();
    155151
    156                   if (!nbsubdomains){
    157                           _error_("nbsubdomains==0");
    158                   }
    159                   if (!subdomains[0].head || !subdomains[0].head->link){
    160                           _error_("!subdomains[0].head || !subdomains[0].head->link");
    161                   }
    162 
     152                  /*Final checks*/
     153                  _assert_(kt==nbt);
     154                  _assert_(nbsubdomains);
     155                  _assert_(subdomains[0].head && subdomains[0].head->link);
    163156          }
    164157        /*}}}*/
    165158        /*FUNCTION Mesh::Mesh(Mesh & Th,Geometry * pGh,Mesh * pBth,long maxnbv_in) COPY{{{*/
    … …  
    17241717                        long i0 = GetId(triangles[it][VerticesOfTriangularEdge[j][0]]);
    17251718                        long i1 = GetId(triangles[it][VerticesOfTriangularEdge[j][1]]);
    17261719                        k = edge4->SortAndFind(i0,i1);
    1727                         if(k>=0){
    1728                                 subdomains[i].direction = (vertices + i0 == edges[k].v[0]) ? 1 : -1;
    1729                                 subdomains[i].edge = edges+k;
    1730                                 Gh.subdomains[i].edge = Gh.edges + k;
    1731                                 Gh.subdomains[i].direction  =  subdomains[i].direction;
    1732                                 Gh.subdomains[i].ReferenceNumber =  subdomains[i].ReferenceNumber;
    1733                         }
    1734                         else
    1735                          _error_("%i should be >=0");
    1736                   }
     1720                        _assert_(k>=0);
     1721                        subdomains[i].direction = (vertices + i0 == edges[k].v[0]) ? 1 : -1;
     1722                        subdomains[i].edge = edges+k;
     1723                        Gh.subdomains[i].edge = Gh.edges + k;
     1724                        Gh.subdomains[i].direction  =  subdomains[i].direction;
     1725                        Gh.subdomains[i].ReferenceNumber =  subdomains[i].ReferenceNumber;
     1726                }
    17371727
    17381728                delete edge4;
    17391729                delete [] colorV;

  • ../trunk-jpl/src/c/modules/InputArtificialNoisex/InputArtificialNoisex.h

     
    1 /*!\file:  InputArtificialNoisex.h
    2  * \brief header file for updating datasets from inputs
    3  */
    4 
    5 #ifndef _INPUTARTIFICIALNOISEX_H
    6 #define _INPUTARTIFICIALNOISEX_H
    7 
    8 #include "../../classes/classes.h"
    9 
    10 /* local prototypes: */
    11 void    InputArtificialNoisex( Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads, Materials* materials,  Parameters* parameters,int enum_name,IssmDouble min,IssmDouble max);
    12 
    13 #endif  /* _UPDATEINPUTSFROMVECTORXX_H */

  • ../trunk-jpl/src/c/modules/InputArtificialNoisex/InputArtificialNoisex.cpp

     
    1 /*!\file InputArtificialNoisex
    2  * \brief: update datasets using  parameter inputs
    3  */
    4 
    5 #include "./InputArtificialNoisex.h"
    6 #include "../../shared/shared.h"
    7 #include "../../toolkits/toolkits.h"
    8 
    9 void InputArtificialNoisex( Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters,int enum_name,IssmDouble min,IssmDouble max){
    10 
    11         int i;
    12 
    13         /*Update elements, nodes, loads and materials from inputs: */
    14         for(i=0;i<elements->Size();i++){
    15                 Element* element=dynamic_cast<Element*>(elements->GetObjectByOffset(i));
    16                 element->InputArtificialNoise(enum_name,min,max);
    17         }
    18 }

  • ../trunk-jpl/src/c/modules/modules.h

     
    5050#include "./InputUpdateFromVectorx/InputUpdateFromVectorx.h"
    5151#include "./InputUpdateFromVectorDakotax/InputUpdateFromVectorDakotax.h"
    5252#include "./InputUpdateFromMatrixDakotax/InputUpdateFromMatrixDakotax.h"
    53 #include "./InputArtificialNoisex/InputArtificialNoisex.h"
    5453#include "./IoModelToConstraintsx/IoModelToConstraintsx.h"
    5554#include "./KMLFileReadx/KMLFileReadx.h"
    5655#include "./KMLMeshWritex/KMLMeshWritex.h"

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

     
    297297                                        ./modules/GetVectorFromInputsx/GetVectorFromInputsx.cpp\
    298298                                        ./modules/InputUpdateFromVectorx/InputUpdateFromVectorx.h\
    299299                                        ./modules/InputUpdateFromVectorx/InputUpdateFromVectorx.cpp\
    300                                         ./modules/InputArtificialNoisex/InputArtificialNoisex.h\
    301                                         ./modules/InputArtificialNoisex/InputArtificialNoisex.cpp\
    302300                                        ./modules/ConfigureObjectsx/ConfigureObjectsx.h\
    303301                                        ./modules/ConfigureObjectsx/ConfigureObjectsx.cpp\
    304302                                        ./modules/ComputeBasalStressx/ComputeBasalStressx.h\

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

     
    5656                          $ENV{ISSM_DIR}/src/c/modules/GetSolutionFromInputsx
    5757                          $ENV{ISSM_DIR}/src/c/modules/GetVectorFromInputsx
    5858                          $ENV{ISSM_DIR}/src/c/modules/InputUpdateFromVectorx
    59                           $ENV{ISSM_DIR}/src/c/modules/InputArtificialNoisex
    6059                          $ENV{ISSM_DIR}/src/c/modules/ConfigureObjectsx
    6160                          $ENV{ISSM_DIR}/src/c/modules/ComputeBasalStressx
    6261                          $ENV{ISSM_DIR}/src/c/modules/ComputeStrainRatex
    … …  
    205204                                        ./modules/GetSolutionFromInputsx/GetSolutionFromInputsx.cpp
    206205                                        ./modules/GetVectorFromInputsx/GetVectorFromInputsx.cpp
    207206                                        ./modules/InputUpdateFromVectorx/InputUpdateFromVectorx.cpp
    208                                         ./modules/InputArtificialNoisex/InputArtificialNoisex.cpp
    209207                                        ./modules/ConfigureObjectsx/ConfigureObjectsx.cpp
    210208                                        ./modules/ComputeBasalStressx/ComputeBasalStressx.cpp
    211209                                        ./modules/ComputeStrainRatex/ComputeStrainRatex.cpp

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

     
    8484                bool  InAnalysis(int analysis_type);
    8585                /*}}}*/
    8686                /*Riftfront specific routines: {{{*/
    87                 bool  PreStable();
    8887                ElementMatrix* PenaltyCreateKMatrixStressbalanceHoriz(IssmDouble kmax);
    8988                ElementVector* PenaltyCreatePVectorStressbalanceHoriz(IssmDouble kmax);
    90                 void  SetPreStable();
    91                 int   PreConstrain(int* punstable);
    9289                int   Constrain(int* punstable);
    9390                void  FreezeConstraints(void);
    9491                bool  IsFrozen(void);
    95                 int   PotentialUnstableConstraint(int* punstable);
    9692                bool  IsInput(int name);
    9793                /*}}}*/
    9894};
    99 
    10095#endif  /* _RIFTFRONT_H_ */

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

     
    732732        else return 0;
    733733}
    734734/*}}}*/
    735 /*FUNCTION Riftfront::PotentialUnstableConstraint {{{*/
    736 int   Riftfront::PotentialUnstableConstraint(int* punstable){
    737 
    738         IssmDouble  penetration;
    739         int         unstable;
    740         IssmDouble  vx1;
    741         IssmDouble  vy1;
    742         IssmDouble  vx2;
    743         IssmDouble  vy2;
    744 
    745         /*Objects: */
    746         Tria       *tria1 = NULL;
    747         Tria       *tria2 = NULL;
    748 
    749         /*enum of element? */
    750         if(elements[0]->ObjectEnum()!=TriaEnum)_error_("only Tria element allowed for Riftfront load!");
    751 
    752         /*recover elements on both side of rift: */
    753         tria1=(Tria*)elements[0];
    754         tria2=(Tria*)elements[1];
    755 
    756         /*First recover velocity: */
    757         tria1->GetInputValue(&vx1,nodes[0],VxEnum);
    758         tria2->GetInputValue(&vx2,nodes[1],VxEnum);
    759         tria1->GetInputValue(&vy1,nodes[0],VyEnum);
    760         tria2->GetInputValue(&vy2,nodes[1],VyEnum);
    761 
    762         /*Node 1 faces node 2, compute penetration of 2 into 1 (V2-V1).N (with N normal vector, and V velocity vector: */
    763         penetration=(vx2-vx1)*normal[0]+(vy2-vy1)*normal[1];
    764 
    765         /*Ok, we are looking for positive penetration in an active constraint: */
    766         if(this->active){
    767                 if (penetration>=0){
    768                         unstable=1;
    769                 }
    770                 else{
    771                         unstable=0;
    772                 }
    773         }
    774         else{
    775                 unstable=0;
    776         }
    777 
    778         /*assign output pointer: */
    779         *punstable=unstable;
    780         return 1;
    781 }
    782 /*}}}*/
    783 /*FUNCTION Riftfront::PreConstrain {{{*/
    784 int   Riftfront::PreConstrain(int* punstable){
    785 
    786         IssmDouble  penetration;
    787         int         unstable;
    788         IssmDouble  vx1;
    789         IssmDouble  vy1;
    790         IssmDouble  vx2;
    791         IssmDouble  vy2;
    792 
    793         /*Objects: */
    794         Tria *tria1 = NULL;
    795         Tria *tria2 = NULL;
    796 
    797         /*enum of element? */
    798         if(elements[0]->ObjectEnum()!=TriaEnum)_error_("only Tria element allowed for Riftfront load!");
    799 
    800         /*recover elements on both side of rift: */
    801         tria1=(Tria*)elements[0];
    802         tria2=(Tria*)elements[1];
    803 
    804         /*First recover velocity: */
    805         tria1->GetInputValue(&vx1,nodes[0],VxEnum);
    806         tria2->GetInputValue(&vx2,nodes[1],VxEnum);
    807         tria1->GetInputValue(&vy1,nodes[0],VyEnum);
    808         tria2->GetInputValue(&vy2,nodes[1],VyEnum);
    809 
    810         /*Node 1 faces node 2, compute penetration of 2 into 1 (V2-V1).N (with N normal vector, and V velocity vector: */
    811         penetration=(vx2-vx1)*normal[0]+(vy2-vy1)*normal[1];
    812 
    813         /*Ok, we are preconstraining here. Ie, anything that penetrates is constrained until stability of the entire set
    814          * of constraints is reached.: */
    815         if(penetration<0){
    816                 if (!this->active){
    817                         /*This is the first time penetration happens: */
    818                         this->active=1;
    819                         unstable=1;
    820                 }
    821                 else{
    822                         /*This constraint was already active: */
    823                         this->active=1;
    824                         unstable=0;
    825                 }
    826         }
    827         else{
    828                 /*No penetration happening. : */
    829                 if (!this->active){
    830                         /*This penalty was not active, and no penetration happening. Do nonthing: */
    831                         this->active=0;
    832                         unstable=0;
    833                 }
    834                 else{
    835                         /*Ok, this penalty wants to get released. But not now, this is preconstraint, not constraint: */
    836                         this->active=1;
    837                         unstable=0;
    838                 }
    839         }
    840 
    841         /*assign output pointer: */
    842         *punstable=unstable;
    843         return 1;
    844 }
    845 /*}}}*/
    846 /*FUNCTION Riftfront::PreStable {{{*/
    847 bool  Riftfront::PreStable(){
    848         return prestable;
    849 }
    850 /*}}}*/
    851 /*FUNCTION Riftfront::SetPreStable {{{*/
    852 void Riftfront::SetPreStable(){
    853         prestable=1;
    854 }
    855 /*}}}*/
    856735/*FUNCTION Riftfront::IsInput{{{*/
    857736bool Riftfront::IsInput(int name){
    858737        if (
Note: See TracBrowser for help on using the repository browser.