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

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

NEW: added Archive/16133-16554
File size: 14.6 KB
RevLine 
[16556]1Index: ../trunk-jpl/src/c/bamg/Geometry.h
2===================================================================
3--- ../trunk-jpl/src/c/bamg/Geometry.h (revision 16236)
4+++ ../trunk-jpl/src/c/bamg/Geometry.h (revision 16237)
5@@ -59,8 +59,7 @@
6 long GetId(const GeomEdge *t) const;
7 long GetId(const Curve *c) const;
8 void UnMarkEdges();
9- GeomEdge *ProjectOnCurve(const Edge &,double,BamgVertex &,VertexOnGeom &) const;
10- GeomEdge *Containing(const R2 P, GeomEdge *start) const;
11+ GeomEdge *ProjectOnCurve(const Edge &,double,BamgVertex &,VertexOnGeom &) const;
12 void WriteGeometry(BamgGeom *bamggeom, BamgOpts*bamgopts);
13 };
14
15Index: ../trunk-jpl/src/c/bamg/Triangle.cpp
16===================================================================
17--- ../trunk-jpl/src/c/bamg/Triangle.cpp (revision 16236)
18+++ ../trunk-jpl/src/c/bamg/Triangle.cpp (revision 16237)
19@@ -251,10 +251,6 @@
20 t->AdjEdgeIndex[aat]=a;
21 }
22 }/*}}}*/
23- /*FUNCTION Triangle::SetAllFlag{{{*/
24- void Triangle::SetAllFlag(int a,int f){
25- AdjEdgeIndex[a] = (AdjEdgeIndex[a] &3) + (1020 & f);
26- }/*}}}*/
27 /*FUNCTION Triangle::SetHidden{{{*/
28 void Triangle::SetHidden(int a){
29 //Get Adjacent Triangle number a
30Index: ../trunk-jpl/src/c/bamg/Geometry.cpp
31===================================================================
32--- ../trunk-jpl/src/c/bamg/Geometry.cpp (revision 16236)
33+++ ../trunk-jpl/src/c/bamg/Geometry.cpp (revision 16237)
34@@ -462,32 +462,6 @@
35 long Geometry::GetId(const Curve * c) const {
36 return c - curves;
37 }/*}}}*/
38- /*FUNCTION Geometry::Containing{{{*/
39- GeomEdge* Geometry::Containing(const R2 P, GeomEdge * start) const {
40- /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshGeom.cpp/Contening)*/
41-
42- GeomEdge* on =start,* pon=0;
43- // walk with the cos on geometry
44- int counter=0;
45- while(pon != on){
46- counter++;
47- _assert_(counter<100);
48- pon = on;
49- R2 A= (*on)[0];
50- R2 B= (*on)[1];
51- R2 AB = B-A;
52- R2 AP = P-A;
53- R2 BP = P-B;
54- if ( (AB,AP) < 0)
55- on = on->Adj[0];
56- else if ( (AB,BP) > 0)
57- on = on->Adj[1];
58- else
59- return on;
60- }
61- return on;
62- }
63- /*}}}*/
64 /*FUNCTION Geometry::PostRead{{{*/
65 void Geometry::PostRead(){
66 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshGeom.cpp/AfterRead)*/
67Index: ../trunk-jpl/src/c/bamg/Triangle.h
68===================================================================
69--- ../trunk-jpl/src/c/bamg/Triangle.h (revision 16236)
70+++ ../trunk-jpl/src/c/bamg/Triangle.h (revision 16237)
71@@ -45,7 +45,6 @@
72 int Locked(int a)const;
73 int Hidden(int a)const;
74 int GetAllflag(int a);
75- void SetAllFlag(int a,int f);
76 double QualityQuad(int a,int option=1) const;
77 short NuEdgeTriangleAdj(int i) const;
78 AdjacentTriangle Adj(int i) const;
79Index: ../trunk-jpl/src/c/bamg/Metric.h
80===================================================================
81--- ../trunk-jpl/src/c/bamg/Metric.h (revision 16236)
82+++ ../trunk-jpl/src/c/bamg/Metric.h (revision 16237)
83@@ -76,8 +76,6 @@
84 double lmax() const;
85 double lmin() const;
86 double Aniso2() const;
87- double Aniso() const;
88- void BoundAniso(const double c);
89 inline void BoundAniso2(const double coef);
90
91 //operators
92Index: ../trunk-jpl/src/c/bamg/EigenMetric.cpp
93===================================================================
94--- ../trunk-jpl/src/c/bamg/EigenMetric.cpp (revision 16236)
95+++ ../trunk-jpl/src/c/bamg/EigenMetric.cpp (revision 16237)
96@@ -98,18 +98,10 @@
97 void EigenMetric::Abs(){
98 lambda1=bamg::Abs(lambda1),lambda2=bamg::Abs(lambda2);
99 }/*}}}*/
100- /*FUNCTION EigenMetric::Aniso{{{*/
101- double EigenMetric::Aniso() const {
102- return sqrt( Aniso2());
103- }/*}}}*/
104 /*FUNCTION EigenMetric::Aniso2{{{*/
105 double EigenMetric::Aniso2() const {
106 return lmax()/lmin();
107 }/*}}}*/
108- /*FUNCTION EigenMetric::BoundAniso{{{*/
109- void EigenMetric::BoundAniso(const double c){
110- BoundAniso2(1/(c*c));
111- }/*}}}*/
112 /*FUNCTION EigenMetric::Echo {{{*/
113 void EigenMetric::Echo(void){
114
115Index: ../trunk-jpl/src/c/bamg/Mesh.cpp
116===================================================================
117--- ../trunk-jpl/src/c/bamg/Mesh.cpp (revision 16236)
118+++ ../trunk-jpl/src/c/bamg/Mesh.cpp (revision 16237)
119@@ -138,28 +138,21 @@
120 triangles[nbt].color = Tho.subdomains[reft[i]].ReferenceNumber;
121 nbt++;
122 }
123- if (kt!=nbt){
124- _error_("kt!=nbt");
125- }
126 if (nbt==0 && nbv==0) {
127 _error_("All triangles have been removed");
128 }
129 delete [] kk;
130 delete [] reft;
131 delete [] refv;
132- //double cutoffradian = 10.0/180.0*Pi;
133 BuildGeometryFromMesh(bamgopts);
134 Gh.PostRead();
135 SetIntCoor();
136 ReconstructExistingMesh();
137
138- if (!nbsubdomains){
139- _error_("nbsubdomains==0");
140- }
141- if (!subdomains[0].head || !subdomains[0].head->link){
142- _error_("!subdomains[0].head || !subdomains[0].head->link");
143- }
144-
145+ /*Final checks*/
146+ _assert_(kt==nbt);
147+ _assert_(nbsubdomains);
148+ _assert_(subdomains[0].head && subdomains[0].head->link);
149 }
150 /*}}}*/
151 /*FUNCTION Mesh::Mesh(Mesh & Th,Geometry * pGh,Mesh * pBth,long maxnbv_in) COPY{{{*/
152@@ -1724,16 +1717,13 @@
153 long i0 = GetId(triangles[it][VerticesOfTriangularEdge[j][0]]);
154 long i1 = GetId(triangles[it][VerticesOfTriangularEdge[j][1]]);
155 k = edge4->SortAndFind(i0,i1);
156- if(k>=0){
157- subdomains[i].direction = (vertices + i0 == edges[k].v[0]) ? 1 : -1;
158- subdomains[i].edge = edges+k;
159- Gh.subdomains[i].edge = Gh.edges + k;
160- Gh.subdomains[i].direction = subdomains[i].direction;
161- Gh.subdomains[i].ReferenceNumber = subdomains[i].ReferenceNumber;
162- }
163- else
164- _error_("%i should be >=0");
165- }
166+ _assert_(k>=0);
167+ subdomains[i].direction = (vertices + i0 == edges[k].v[0]) ? 1 : -1;
168+ subdomains[i].edge = edges+k;
169+ Gh.subdomains[i].edge = Gh.edges + k;
170+ Gh.subdomains[i].direction = subdomains[i].direction;
171+ Gh.subdomains[i].ReferenceNumber = subdomains[i].ReferenceNumber;
172+ }
173
174 delete edge4;
175 delete [] colorV;
176Index: ../trunk-jpl/src/c/modules/InputArtificialNoisex/InputArtificialNoisex.h
177===================================================================
178--- ../trunk-jpl/src/c/modules/InputArtificialNoisex/InputArtificialNoisex.h (revision 16236)
179+++ ../trunk-jpl/src/c/modules/InputArtificialNoisex/InputArtificialNoisex.h (revision 16237)
180@@ -1,13 +0,0 @@
181-/*!\file: InputArtificialNoisex.h
182- * \brief header file for updating datasets from inputs
183- */
184-
185-#ifndef _INPUTARTIFICIALNOISEX_H
186-#define _INPUTARTIFICIALNOISEX_H
187-
188-#include "../../classes/classes.h"
189-
190-/* local prototypes: */
191-void InputArtificialNoisex( Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads, Materials* materials, Parameters* parameters,int enum_name,IssmDouble min,IssmDouble max);
192-
193-#endif /* _UPDATEINPUTSFROMVECTORXX_H */
194Index: ../trunk-jpl/src/c/modules/InputArtificialNoisex/InputArtificialNoisex.cpp
195===================================================================
196--- ../trunk-jpl/src/c/modules/InputArtificialNoisex/InputArtificialNoisex.cpp (revision 16236)
197+++ ../trunk-jpl/src/c/modules/InputArtificialNoisex/InputArtificialNoisex.cpp (revision 16237)
198@@ -1,18 +0,0 @@
199-/*!\file InputArtificialNoisex
200- * \brief: update datasets using parameter inputs
201- */
202-
203-#include "./InputArtificialNoisex.h"
204-#include "../../shared/shared.h"
205-#include "../../toolkits/toolkits.h"
206-
207-void InputArtificialNoisex( Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters,int enum_name,IssmDouble min,IssmDouble max){
208-
209- int i;
210-
211- /*Update elements, nodes, loads and materials from inputs: */
212- for(i=0;i<elements->Size();i++){
213- Element* element=dynamic_cast<Element*>(elements->GetObjectByOffset(i));
214- element->InputArtificialNoise(enum_name,min,max);
215- }
216-}
217Index: ../trunk-jpl/src/c/modules/modules.h
218===================================================================
219--- ../trunk-jpl/src/c/modules/modules.h (revision 16236)
220+++ ../trunk-jpl/src/c/modules/modules.h (revision 16237)
221@@ -50,7 +50,6 @@
222 #include "./InputUpdateFromVectorx/InputUpdateFromVectorx.h"
223 #include "./InputUpdateFromVectorDakotax/InputUpdateFromVectorDakotax.h"
224 #include "./InputUpdateFromMatrixDakotax/InputUpdateFromMatrixDakotax.h"
225-#include "./InputArtificialNoisex/InputArtificialNoisex.h"
226 #include "./IoModelToConstraintsx/IoModelToConstraintsx.h"
227 #include "./KMLFileReadx/KMLFileReadx.h"
228 #include "./KMLMeshWritex/KMLMeshWritex.h"
229Index: ../trunk-jpl/src/c/Makefile.am
230===================================================================
231--- ../trunk-jpl/src/c/Makefile.am (revision 16236)
232+++ ../trunk-jpl/src/c/Makefile.am (revision 16237)
233@@ -297,8 +297,6 @@
234 ./modules/GetVectorFromInputsx/GetVectorFromInputsx.cpp\
235 ./modules/InputUpdateFromVectorx/InputUpdateFromVectorx.h\
236 ./modules/InputUpdateFromVectorx/InputUpdateFromVectorx.cpp\
237- ./modules/InputArtificialNoisex/InputArtificialNoisex.h\
238- ./modules/InputArtificialNoisex/InputArtificialNoisex.cpp\
239 ./modules/ConfigureObjectsx/ConfigureObjectsx.h\
240 ./modules/ConfigureObjectsx/ConfigureObjectsx.cpp\
241 ./modules/ComputeBasalStressx/ComputeBasalStressx.h\
242Index: ../trunk-jpl/src/c/CMakeLists.txt
243===================================================================
244--- ../trunk-jpl/src/c/CMakeLists.txt (revision 16236)
245+++ ../trunk-jpl/src/c/CMakeLists.txt (revision 16237)
246@@ -56,7 +56,6 @@
247 $ENV{ISSM_DIR}/src/c/modules/GetSolutionFromInputsx
248 $ENV{ISSM_DIR}/src/c/modules/GetVectorFromInputsx
249 $ENV{ISSM_DIR}/src/c/modules/InputUpdateFromVectorx
250- $ENV{ISSM_DIR}/src/c/modules/InputArtificialNoisex
251 $ENV{ISSM_DIR}/src/c/modules/ConfigureObjectsx
252 $ENV{ISSM_DIR}/src/c/modules/ComputeBasalStressx
253 $ENV{ISSM_DIR}/src/c/modules/ComputeStrainRatex
254@@ -205,7 +204,6 @@
255 ./modules/GetSolutionFromInputsx/GetSolutionFromInputsx.cpp
256 ./modules/GetVectorFromInputsx/GetVectorFromInputsx.cpp
257 ./modules/InputUpdateFromVectorx/InputUpdateFromVectorx.cpp
258- ./modules/InputArtificialNoisex/InputArtificialNoisex.cpp
259 ./modules/ConfigureObjectsx/ConfigureObjectsx.cpp
260 ./modules/ComputeBasalStressx/ComputeBasalStressx.cpp
261 ./modules/ComputeStrainRatex/ComputeStrainRatex.cpp
262Index: ../trunk-jpl/src/c/classes/Loads/Riftfront.h
263===================================================================
264--- ../trunk-jpl/src/c/classes/Loads/Riftfront.h (revision 16236)
265+++ ../trunk-jpl/src/c/classes/Loads/Riftfront.h (revision 16237)
266@@ -84,17 +84,12 @@
267 bool InAnalysis(int analysis_type);
268 /*}}}*/
269 /*Riftfront specific routines: {{{*/
270- bool PreStable();
271 ElementMatrix* PenaltyCreateKMatrixStressbalanceHoriz(IssmDouble kmax);
272 ElementVector* PenaltyCreatePVectorStressbalanceHoriz(IssmDouble kmax);
273- void SetPreStable();
274- int PreConstrain(int* punstable);
275 int Constrain(int* punstable);
276 void FreezeConstraints(void);
277 bool IsFrozen(void);
278- int PotentialUnstableConstraint(int* punstable);
279 bool IsInput(int name);
280 /*}}}*/
281 };
282-
283 #endif /* _RIFTFRONT_H_ */
284Index: ../trunk-jpl/src/c/classes/Loads/Riftfront.cpp
285===================================================================
286--- ../trunk-jpl/src/c/classes/Loads/Riftfront.cpp (revision 16236)
287+++ ../trunk-jpl/src/c/classes/Loads/Riftfront.cpp (revision 16237)
288@@ -732,127 +732,6 @@
289 else return 0;
290 }
291 /*}}}*/
292-/*FUNCTION Riftfront::PotentialUnstableConstraint {{{*/
293-int Riftfront::PotentialUnstableConstraint(int* punstable){
294-
295- IssmDouble penetration;
296- int unstable;
297- IssmDouble vx1;
298- IssmDouble vy1;
299- IssmDouble vx2;
300- IssmDouble vy2;
301-
302- /*Objects: */
303- Tria *tria1 = NULL;
304- Tria *tria2 = NULL;
305-
306- /*enum of element? */
307- if(elements[0]->ObjectEnum()!=TriaEnum)_error_("only Tria element allowed for Riftfront load!");
308-
309- /*recover elements on both side of rift: */
310- tria1=(Tria*)elements[0];
311- tria2=(Tria*)elements[1];
312-
313- /*First recover velocity: */
314- tria1->GetInputValue(&vx1,nodes[0],VxEnum);
315- tria2->GetInputValue(&vx2,nodes[1],VxEnum);
316- tria1->GetInputValue(&vy1,nodes[0],VyEnum);
317- tria2->GetInputValue(&vy2,nodes[1],VyEnum);
318-
319- /*Node 1 faces node 2, compute penetration of 2 into 1 (V2-V1).N (with N normal vector, and V velocity vector: */
320- penetration=(vx2-vx1)*normal[0]+(vy2-vy1)*normal[1];
321-
322- /*Ok, we are looking for positive penetration in an active constraint: */
323- if(this->active){
324- if (penetration>=0){
325- unstable=1;
326- }
327- else{
328- unstable=0;
329- }
330- }
331- else{
332- unstable=0;
333- }
334-
335- /*assign output pointer: */
336- *punstable=unstable;
337- return 1;
338-}
339-/*}}}*/
340-/*FUNCTION Riftfront::PreConstrain {{{*/
341-int Riftfront::PreConstrain(int* punstable){
342-
343- IssmDouble penetration;
344- int unstable;
345- IssmDouble vx1;
346- IssmDouble vy1;
347- IssmDouble vx2;
348- IssmDouble vy2;
349-
350- /*Objects: */
351- Tria *tria1 = NULL;
352- Tria *tria2 = NULL;
353-
354- /*enum of element? */
355- if(elements[0]->ObjectEnum()!=TriaEnum)_error_("only Tria element allowed for Riftfront load!");
356-
357- /*recover elements on both side of rift: */
358- tria1=(Tria*)elements[0];
359- tria2=(Tria*)elements[1];
360-
361- /*First recover velocity: */
362- tria1->GetInputValue(&vx1,nodes[0],VxEnum);
363- tria2->GetInputValue(&vx2,nodes[1],VxEnum);
364- tria1->GetInputValue(&vy1,nodes[0],VyEnum);
365- tria2->GetInputValue(&vy2,nodes[1],VyEnum);
366-
367- /*Node 1 faces node 2, compute penetration of 2 into 1 (V2-V1).N (with N normal vector, and V velocity vector: */
368- penetration=(vx2-vx1)*normal[0]+(vy2-vy1)*normal[1];
369-
370- /*Ok, we are preconstraining here. Ie, anything that penetrates is constrained until stability of the entire set
371- * of constraints is reached.: */
372- if(penetration<0){
373- if (!this->active){
374- /*This is the first time penetration happens: */
375- this->active=1;
376- unstable=1;
377- }
378- else{
379- /*This constraint was already active: */
380- this->active=1;
381- unstable=0;
382- }
383- }
384- else{
385- /*No penetration happening. : */
386- if (!this->active){
387- /*This penalty was not active, and no penetration happening. Do nonthing: */
388- this->active=0;
389- unstable=0;
390- }
391- else{
392- /*Ok, this penalty wants to get released. But not now, this is preconstraint, not constraint: */
393- this->active=1;
394- unstable=0;
395- }
396- }
397-
398- /*assign output pointer: */
399- *punstable=unstable;
400- return 1;
401-}
402-/*}}}*/
403-/*FUNCTION Riftfront::PreStable {{{*/
404-bool Riftfront::PreStable(){
405- return prestable;
406-}
407-/*}}}*/
408-/*FUNCTION Riftfront::SetPreStable {{{*/
409-void Riftfront::SetPreStable(){
410- prestable=1;
411-}
412-/*}}}*/
413 /*FUNCTION Riftfront::IsInput{{{*/
414 bool Riftfront::IsInput(int name){
415 if (
Note: See TracBrowser for help on using the repository browser.