Index: ../trunk-jpl/src/c/bamg/Geometry.h =================================================================== --- ../trunk-jpl/src/c/bamg/Geometry.h (revision 16236) +++ ../trunk-jpl/src/c/bamg/Geometry.h (revision 16237) @@ -59,8 +59,7 @@ long GetId(const GeomEdge *t) const; long GetId(const Curve *c) const; void UnMarkEdges(); - GeomEdge *ProjectOnCurve(const Edge &,double,BamgVertex &,VertexOnGeom &) const; - GeomEdge *Containing(const R2 P, GeomEdge *start) const; + GeomEdge *ProjectOnCurve(const Edge &,double,BamgVertex &,VertexOnGeom &) const; void WriteGeometry(BamgGeom *bamggeom, BamgOpts*bamgopts); }; Index: ../trunk-jpl/src/c/bamg/Triangle.cpp =================================================================== --- ../trunk-jpl/src/c/bamg/Triangle.cpp (revision 16236) +++ ../trunk-jpl/src/c/bamg/Triangle.cpp (revision 16237) @@ -251,10 +251,6 @@ t->AdjEdgeIndex[aat]=a; } }/*}}}*/ - /*FUNCTION Triangle::SetAllFlag{{{*/ - void Triangle::SetAllFlag(int a,int f){ - AdjEdgeIndex[a] = (AdjEdgeIndex[a] &3) + (1020 & f); - }/*}}}*/ /*FUNCTION Triangle::SetHidden{{{*/ void Triangle::SetHidden(int a){ //Get Adjacent Triangle number a Index: ../trunk-jpl/src/c/bamg/Geometry.cpp =================================================================== --- ../trunk-jpl/src/c/bamg/Geometry.cpp (revision 16236) +++ ../trunk-jpl/src/c/bamg/Geometry.cpp (revision 16237) @@ -462,32 +462,6 @@ long Geometry::GetId(const Curve * c) const { return c - curves; }/*}}}*/ - /*FUNCTION Geometry::Containing{{{*/ - GeomEdge* Geometry::Containing(const R2 P, GeomEdge * start) const { - /*Original code from Frederic Hecht (BAMG v1.01, MeshGeom.cpp/Contening)*/ - - GeomEdge* on =start,* pon=0; - // walk with the cos on geometry - int counter=0; - while(pon != on){ - counter++; - _assert_(counter<100); - pon = on; - R2 A= (*on)[0]; - R2 B= (*on)[1]; - R2 AB = B-A; - R2 AP = P-A; - R2 BP = P-B; - if ( (AB,AP) < 0) - on = on->Adj[0]; - else if ( (AB,BP) > 0) - on = on->Adj[1]; - else - return on; - } - return on; - } - /*}}}*/ /*FUNCTION Geometry::PostRead{{{*/ void Geometry::PostRead(){ /*Original code from Frederic Hecht (BAMG v1.01, MeshGeom.cpp/AfterRead)*/ Index: ../trunk-jpl/src/c/bamg/Triangle.h =================================================================== --- ../trunk-jpl/src/c/bamg/Triangle.h (revision 16236) +++ ../trunk-jpl/src/c/bamg/Triangle.h (revision 16237) @@ -45,7 +45,6 @@ int Locked(int a)const; int Hidden(int a)const; int GetAllflag(int a); - void SetAllFlag(int a,int f); double QualityQuad(int a,int option=1) const; short NuEdgeTriangleAdj(int i) const; AdjacentTriangle Adj(int i) const; Index: ../trunk-jpl/src/c/bamg/Metric.h =================================================================== --- ../trunk-jpl/src/c/bamg/Metric.h (revision 16236) +++ ../trunk-jpl/src/c/bamg/Metric.h (revision 16237) @@ -76,8 +76,6 @@ double lmax() const; double lmin() const; double Aniso2() const; - double Aniso() const; - void BoundAniso(const double c); inline void BoundAniso2(const double coef); //operators Index: ../trunk-jpl/src/c/bamg/EigenMetric.cpp =================================================================== --- ../trunk-jpl/src/c/bamg/EigenMetric.cpp (revision 16236) +++ ../trunk-jpl/src/c/bamg/EigenMetric.cpp (revision 16237) @@ -98,18 +98,10 @@ void EigenMetric::Abs(){ lambda1=bamg::Abs(lambda1),lambda2=bamg::Abs(lambda2); }/*}}}*/ - /*FUNCTION EigenMetric::Aniso{{{*/ - double EigenMetric::Aniso() const { - return sqrt( Aniso2()); - }/*}}}*/ /*FUNCTION EigenMetric::Aniso2{{{*/ double EigenMetric::Aniso2() const { return lmax()/lmin(); }/*}}}*/ - /*FUNCTION EigenMetric::BoundAniso{{{*/ - void EigenMetric::BoundAniso(const double c){ - BoundAniso2(1/(c*c)); - }/*}}}*/ /*FUNCTION EigenMetric::Echo {{{*/ void EigenMetric::Echo(void){ Index: ../trunk-jpl/src/c/bamg/Mesh.cpp =================================================================== --- ../trunk-jpl/src/c/bamg/Mesh.cpp (revision 16236) +++ ../trunk-jpl/src/c/bamg/Mesh.cpp (revision 16237) @@ -138,28 +138,21 @@ triangles[nbt].color = Tho.subdomains[reft[i]].ReferenceNumber; nbt++; } - if (kt!=nbt){ - _error_("kt!=nbt"); - } if (nbt==0 && nbv==0) { _error_("All triangles have been removed"); } delete [] kk; delete [] reft; delete [] refv; - //double cutoffradian = 10.0/180.0*Pi; BuildGeometryFromMesh(bamgopts); Gh.PostRead(); SetIntCoor(); ReconstructExistingMesh(); - if (!nbsubdomains){ - _error_("nbsubdomains==0"); - } - if (!subdomains[0].head || !subdomains[0].head->link){ - _error_("!subdomains[0].head || !subdomains[0].head->link"); - } - + /*Final checks*/ + _assert_(kt==nbt); + _assert_(nbsubdomains); + _assert_(subdomains[0].head && subdomains[0].head->link); } /*}}}*/ /*FUNCTION Mesh::Mesh(Mesh & Th,Geometry * pGh,Mesh * pBth,long maxnbv_in) COPY{{{*/ @@ -1724,16 +1717,13 @@ long i0 = GetId(triangles[it][VerticesOfTriangularEdge[j][0]]); long i1 = GetId(triangles[it][VerticesOfTriangularEdge[j][1]]); k = edge4->SortAndFind(i0,i1); - if(k>=0){ - subdomains[i].direction = (vertices + i0 == edges[k].v[0]) ? 1 : -1; - subdomains[i].edge = edges+k; - Gh.subdomains[i].edge = Gh.edges + k; - Gh.subdomains[i].direction = subdomains[i].direction; - Gh.subdomains[i].ReferenceNumber = subdomains[i].ReferenceNumber; - } - else - _error_("%i should be >=0"); - } + _assert_(k>=0); + subdomains[i].direction = (vertices + i0 == edges[k].v[0]) ? 1 : -1; + subdomains[i].edge = edges+k; + Gh.subdomains[i].edge = Gh.edges + k; + Gh.subdomains[i].direction = subdomains[i].direction; + Gh.subdomains[i].ReferenceNumber = subdomains[i].ReferenceNumber; + } delete edge4; delete [] colorV; Index: ../trunk-jpl/src/c/modules/InputArtificialNoisex/InputArtificialNoisex.h =================================================================== --- ../trunk-jpl/src/c/modules/InputArtificialNoisex/InputArtificialNoisex.h (revision 16236) +++ ../trunk-jpl/src/c/modules/InputArtificialNoisex/InputArtificialNoisex.h (revision 16237) @@ -1,13 +0,0 @@ -/*!\file: InputArtificialNoisex.h - * \brief header file for updating datasets from inputs - */ - -#ifndef _INPUTARTIFICIALNOISEX_H -#define _INPUTARTIFICIALNOISEX_H - -#include "../../classes/classes.h" - -/* local prototypes: */ -void InputArtificialNoisex( Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads, Materials* materials, Parameters* parameters,int enum_name,IssmDouble min,IssmDouble max); - -#endif /* _UPDATEINPUTSFROMVECTORXX_H */ Index: ../trunk-jpl/src/c/modules/InputArtificialNoisex/InputArtificialNoisex.cpp =================================================================== --- ../trunk-jpl/src/c/modules/InputArtificialNoisex/InputArtificialNoisex.cpp (revision 16236) +++ ../trunk-jpl/src/c/modules/InputArtificialNoisex/InputArtificialNoisex.cpp (revision 16237) @@ -1,18 +0,0 @@ -/*!\file InputArtificialNoisex - * \brief: update datasets using parameter inputs - */ - -#include "./InputArtificialNoisex.h" -#include "../../shared/shared.h" -#include "../../toolkits/toolkits.h" - -void InputArtificialNoisex( Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters,int enum_name,IssmDouble min,IssmDouble max){ - - int i; - - /*Update elements, nodes, loads and materials from inputs: */ - for(i=0;iSize();i++){ - Element* element=dynamic_cast(elements->GetObjectByOffset(i)); - element->InputArtificialNoise(enum_name,min,max); - } -} Index: ../trunk-jpl/src/c/modules/modules.h =================================================================== --- ../trunk-jpl/src/c/modules/modules.h (revision 16236) +++ ../trunk-jpl/src/c/modules/modules.h (revision 16237) @@ -50,7 +50,6 @@ #include "./InputUpdateFromVectorx/InputUpdateFromVectorx.h" #include "./InputUpdateFromVectorDakotax/InputUpdateFromVectorDakotax.h" #include "./InputUpdateFromMatrixDakotax/InputUpdateFromMatrixDakotax.h" -#include "./InputArtificialNoisex/InputArtificialNoisex.h" #include "./IoModelToConstraintsx/IoModelToConstraintsx.h" #include "./KMLFileReadx/KMLFileReadx.h" #include "./KMLMeshWritex/KMLMeshWritex.h" Index: ../trunk-jpl/src/c/Makefile.am =================================================================== --- ../trunk-jpl/src/c/Makefile.am (revision 16236) +++ ../trunk-jpl/src/c/Makefile.am (revision 16237) @@ -297,8 +297,6 @@ ./modules/GetVectorFromInputsx/GetVectorFromInputsx.cpp\ ./modules/InputUpdateFromVectorx/InputUpdateFromVectorx.h\ ./modules/InputUpdateFromVectorx/InputUpdateFromVectorx.cpp\ - ./modules/InputArtificialNoisex/InputArtificialNoisex.h\ - ./modules/InputArtificialNoisex/InputArtificialNoisex.cpp\ ./modules/ConfigureObjectsx/ConfigureObjectsx.h\ ./modules/ConfigureObjectsx/ConfigureObjectsx.cpp\ ./modules/ComputeBasalStressx/ComputeBasalStressx.h\ Index: ../trunk-jpl/src/c/CMakeLists.txt =================================================================== --- ../trunk-jpl/src/c/CMakeLists.txt (revision 16236) +++ ../trunk-jpl/src/c/CMakeLists.txt (revision 16237) @@ -56,7 +56,6 @@ $ENV{ISSM_DIR}/src/c/modules/GetSolutionFromInputsx $ENV{ISSM_DIR}/src/c/modules/GetVectorFromInputsx $ENV{ISSM_DIR}/src/c/modules/InputUpdateFromVectorx - $ENV{ISSM_DIR}/src/c/modules/InputArtificialNoisex $ENV{ISSM_DIR}/src/c/modules/ConfigureObjectsx $ENV{ISSM_DIR}/src/c/modules/ComputeBasalStressx $ENV{ISSM_DIR}/src/c/modules/ComputeStrainRatex @@ -205,7 +204,6 @@ ./modules/GetSolutionFromInputsx/GetSolutionFromInputsx.cpp ./modules/GetVectorFromInputsx/GetVectorFromInputsx.cpp ./modules/InputUpdateFromVectorx/InputUpdateFromVectorx.cpp - ./modules/InputArtificialNoisex/InputArtificialNoisex.cpp ./modules/ConfigureObjectsx/ConfigureObjectsx.cpp ./modules/ComputeBasalStressx/ComputeBasalStressx.cpp ./modules/ComputeStrainRatex/ComputeStrainRatex.cpp Index: ../trunk-jpl/src/c/classes/Loads/Riftfront.h =================================================================== --- ../trunk-jpl/src/c/classes/Loads/Riftfront.h (revision 16236) +++ ../trunk-jpl/src/c/classes/Loads/Riftfront.h (revision 16237) @@ -84,17 +84,12 @@ bool InAnalysis(int analysis_type); /*}}}*/ /*Riftfront specific routines: {{{*/ - bool PreStable(); ElementMatrix* PenaltyCreateKMatrixStressbalanceHoriz(IssmDouble kmax); ElementVector* PenaltyCreatePVectorStressbalanceHoriz(IssmDouble kmax); - void SetPreStable(); - int PreConstrain(int* punstable); int Constrain(int* punstable); void FreezeConstraints(void); bool IsFrozen(void); - int PotentialUnstableConstraint(int* punstable); bool IsInput(int name); /*}}}*/ }; - #endif /* _RIFTFRONT_H_ */ Index: ../trunk-jpl/src/c/classes/Loads/Riftfront.cpp =================================================================== --- ../trunk-jpl/src/c/classes/Loads/Riftfront.cpp (revision 16236) +++ ../trunk-jpl/src/c/classes/Loads/Riftfront.cpp (revision 16237) @@ -732,127 +732,6 @@ else return 0; } /*}}}*/ -/*FUNCTION Riftfront::PotentialUnstableConstraint {{{*/ -int Riftfront::PotentialUnstableConstraint(int* punstable){ - - IssmDouble penetration; - int unstable; - IssmDouble vx1; - IssmDouble vy1; - IssmDouble vx2; - IssmDouble vy2; - - /*Objects: */ - Tria *tria1 = NULL; - Tria *tria2 = NULL; - - /*enum of element? */ - if(elements[0]->ObjectEnum()!=TriaEnum)_error_("only Tria element allowed for Riftfront load!"); - - /*recover elements on both side of rift: */ - tria1=(Tria*)elements[0]; - tria2=(Tria*)elements[1]; - - /*First recover velocity: */ - tria1->GetInputValue(&vx1,nodes[0],VxEnum); - tria2->GetInputValue(&vx2,nodes[1],VxEnum); - tria1->GetInputValue(&vy1,nodes[0],VyEnum); - tria2->GetInputValue(&vy2,nodes[1],VyEnum); - - /*Node 1 faces node 2, compute penetration of 2 into 1 (V2-V1).N (with N normal vector, and V velocity vector: */ - penetration=(vx2-vx1)*normal[0]+(vy2-vy1)*normal[1]; - - /*Ok, we are looking for positive penetration in an active constraint: */ - if(this->active){ - if (penetration>=0){ - unstable=1; - } - else{ - unstable=0; - } - } - else{ - unstable=0; - } - - /*assign output pointer: */ - *punstable=unstable; - return 1; -} -/*}}}*/ -/*FUNCTION Riftfront::PreConstrain {{{*/ -int Riftfront::PreConstrain(int* punstable){ - - IssmDouble penetration; - int unstable; - IssmDouble vx1; - IssmDouble vy1; - IssmDouble vx2; - IssmDouble vy2; - - /*Objects: */ - Tria *tria1 = NULL; - Tria *tria2 = NULL; - - /*enum of element? */ - if(elements[0]->ObjectEnum()!=TriaEnum)_error_("only Tria element allowed for Riftfront load!"); - - /*recover elements on both side of rift: */ - tria1=(Tria*)elements[0]; - tria2=(Tria*)elements[1]; - - /*First recover velocity: */ - tria1->GetInputValue(&vx1,nodes[0],VxEnum); - tria2->GetInputValue(&vx2,nodes[1],VxEnum); - tria1->GetInputValue(&vy1,nodes[0],VyEnum); - tria2->GetInputValue(&vy2,nodes[1],VyEnum); - - /*Node 1 faces node 2, compute penetration of 2 into 1 (V2-V1).N (with N normal vector, and V velocity vector: */ - penetration=(vx2-vx1)*normal[0]+(vy2-vy1)*normal[1]; - - /*Ok, we are preconstraining here. Ie, anything that penetrates is constrained until stability of the entire set - * of constraints is reached.: */ - if(penetration<0){ - if (!this->active){ - /*This is the first time penetration happens: */ - this->active=1; - unstable=1; - } - else{ - /*This constraint was already active: */ - this->active=1; - unstable=0; - } - } - else{ - /*No penetration happening. : */ - if (!this->active){ - /*This penalty was not active, and no penetration happening. Do nonthing: */ - this->active=0; - unstable=0; - } - else{ - /*Ok, this penalty wants to get released. But not now, this is preconstraint, not constraint: */ - this->active=1; - unstable=0; - } - } - - /*assign output pointer: */ - *punstable=unstable; - return 1; -} -/*}}}*/ -/*FUNCTION Riftfront::PreStable {{{*/ -bool Riftfront::PreStable(){ - return prestable; -} -/*}}}*/ -/*FUNCTION Riftfront::SetPreStable {{{*/ -void Riftfront::SetPreStable(){ - prestable=1; -} -/*}}}*/ /*FUNCTION Riftfront::IsInput{{{*/ bool Riftfront::IsInput(int name){ if (