source: issm/oecreview/Archive/14312-15392/ISSM-14778-14779.diff@ 15393

Last change on this file since 15393 was 15393, checked in by Mathieu Morlighem, 12 years ago

NEW: adding Archive/14312-15392 for oecreview
File size: 5.4 KB
RevLine 
[15393]1Index: ../trunk-jpl/src/c/classes/gauss/GaussTria.cpp
2===================================================================
3--- ../trunk-jpl/src/c/classes/gauss/GaussTria.cpp (revision 14778)
4+++ ../trunk-jpl/src/c/classes/gauss/GaussTria.cpp (revision 14779)
5@@ -92,6 +92,162 @@
6 xDelete<double>(seg_weights);
7 }
8 /*}}}*/
9+/*FUNCTION GaussTria::GaussTria(int index,double r1,double r2,int order) {{{*/
10+GaussTria::GaussTria(int index,IssmDouble r1,IssmDouble r2,bool mainlyfloating,int order){
11+
12+ /*
13+ * ^
14+ * |
15+ * 1|\
16+ * | \
17+ * | \
18+ * | \
19+ * | \
20+ * | \
21+ * | +(x,y) \
22+ * | \
23+ * +---------------+-->
24+ * 0 1
25+ *
26+ */
27+ int ig;
28+ IssmDouble x,y;
29+ IssmDouble xy_list[3][2];
30+
31+ if(mainlyfloating){
32+ /*Get gauss points*/
33+ GaussLegendreTria(&this->numgauss,&this->coords1,&this->coords2,&this->coords3,&this->weights,order);
34+
35+ xy_list[0][0]=0; xy_list[0][1]=0;
36+ xy_list[1][0]=r1; xy_list[1][1]=0;
37+ xy_list[2][0]=0; xy_list[2][1]=r2;
38+
39+ for(ig=0;ig<this->numgauss;ig++){
40+ x = this->coords1[ig]*xy_list[0][0] + this->coords2[ig]*xy_list[1][0] + this->coords3[ig]*xy_list[2][0];
41+ y = this->coords1[ig]*xy_list[0][1] + this->coords2[ig]*xy_list[1][1] + this->coords3[ig]*xy_list[2][1];
42+
43+ switch(index){
44+ case 0:
45+ this->coords1[ig] = 1.-x-y;
46+ this->coords2[ig] = x;
47+ this->coords3[ig] = y;
48+ break;
49+ case 1:
50+ this->coords1[ig] = y;
51+ this->coords2[ig] = 1.-x-y;
52+ this->coords3[ig] = x;
53+ break;
54+ case 2:
55+ this->coords1[ig] = x;
56+ this->coords2[ig] = y;
57+ this->coords3[ig] = 1.-x-y;
58+ break;
59+ default:
60+ _error_("index "<<index<<" not supported yet");
61+ }
62+ this->weights[ig] = this->weights[ig]*r1*r2;
63+ }
64+ }
65+ else{
66+ /*Double number of gauss points*/
67+ GaussTria *gauss1 = NULL;
68+ GaussTria *gauss2 = NULL;
69+ gauss1=new GaussTria(order);
70+ gauss2=new GaussTria(order);
71+
72+ xy_list[0][0]=r1; xy_list[0][1]=0;
73+ xy_list[1][0]=0; xy_list[1][1]=1.;
74+ xy_list[2][0]=0; xy_list[2][1]=r2;
75+
76+ //gauss1->Echo();
77+ for(ig=0;ig<gauss1->numgauss;ig++){
78+ x = gauss1->coords1[ig]*xy_list[0][0] + gauss1->coords2[ig]*xy_list[1][0] + gauss1->coords3[ig]*xy_list[2][0];
79+ y = gauss1->coords1[ig]*xy_list[0][1] + gauss1->coords2[ig]*xy_list[1][1] + gauss1->coords3[ig]*xy_list[2][1];
80+
81+ switch(index){
82+ case 0:
83+ gauss1->coords1[ig] = 1.-x-y;
84+ gauss1->coords2[ig] = x;
85+ gauss1->coords3[ig] = y;
86+ break;
87+ case 1:
88+ gauss1->coords1[ig] = y;
89+ gauss1->coords2[ig] = 1.-x-y;
90+ gauss1->coords3[ig] = x;
91+ break;
92+ case 2:
93+ gauss1->coords1[ig] = x;
94+ gauss1->coords2[ig] = y;
95+ gauss1->coords3[ig] = 1.-x-y;
96+ break;
97+ default:
98+ _error_("index "<<index<<" not supported yet");
99+ }
100+ gauss1->weights[ig] = gauss1->weights[ig]*r1*(1-r2);
101+ }
102+ //gauss1->Echo();
103+ xy_list[0][0]=r1; xy_list[0][1]=0;
104+ xy_list[1][0]=1.; xy_list[1][1]=0;
105+ xy_list[2][0]=0; xy_list[2][1]=1.;
106+
107+ //gauss2->Echo();
108+ for(ig=0;ig<gauss2->numgauss;ig++){
109+ x = gauss2->coords1[ig]*xy_list[0][0] + gauss2->coords2[ig]*xy_list[1][0] + gauss2->coords3[ig]*xy_list[2][0];
110+ y = gauss2->coords1[ig]*xy_list[0][1] + gauss2->coords2[ig]*xy_list[1][1] + gauss2->coords3[ig]*xy_list[2][1];
111+
112+ switch(index){
113+ case 0:
114+ gauss2->coords1[ig] = 1.-x-y;
115+ gauss2->coords2[ig] = x;
116+ gauss2->coords3[ig] = y;
117+ break;
118+ case 1:
119+ gauss2->coords1[ig] = y;
120+ gauss2->coords2[ig] = 1.-x-y;
121+ gauss2->coords3[ig] = x;
122+ break;
123+ case 2:
124+ gauss2->coords1[ig] = x;
125+ gauss2->coords2[ig] = y;
126+ gauss2->coords3[ig] = 1.-x-y;
127+ break;
128+ default:
129+ _error_("index "<<index<<" not supported yet");
130+ }
131+ gauss2->weights[ig] = gauss2->weights[ig]*(1-r1);
132+ }
133+
134+ this->numgauss = gauss1->numgauss + gauss2->numgauss;
135+ this->coords1=xNew<IssmPDouble>(this->numgauss);
136+ this->coords2=xNew<IssmPDouble>(this->numgauss);
137+ this->coords3=xNew<IssmPDouble>(this->numgauss);
138+ this->weights=xNew<IssmPDouble>(this->numgauss);
139+
140+ for(ig=0;ig<gauss1->numgauss;ig++){ // Add the first triangle gauss points
141+ this->coords1[ig]=gauss1->coords1[ig];
142+ this->coords2[ig]=gauss1->coords2[ig];
143+ this->coords3[ig]=gauss1->coords3[ig];
144+ this->weights[ig]=gauss1->weights[ig];
145+ }
146+ for(ig=0;ig<gauss2->numgauss;ig++){ // Add the second triangle gauss points
147+ this->coords1[gauss1->numgauss+ig]=gauss2->coords1[ig];
148+ this->coords2[gauss1->numgauss+ig]=gauss2->coords2[ig];
149+ this->coords3[gauss1->numgauss+ig]=gauss2->coords3[ig];
150+ this->weights[gauss1->numgauss+ig]=gauss2->weights[ig];
151+ }
152+
153+ /*Delete gauss points*/
154+ delete gauss1;
155+ delete gauss2;
156+ }
157+
158+ /*Initialize static fields as undefined*/
159+ weight=UNDEF;
160+ coord1=UNDEF;
161+ coord2=UNDEF;
162+ coord3=UNDEF;
163+}
164+/*}}}*/
165 /*FUNCTION GaussTria::~GaussTria(){{{*/
166 GaussTria::~GaussTria(){
167 xDelete<IssmPDouble>(weights);
168Index: ../trunk-jpl/src/c/classes/gauss/GaussTria.h
169===================================================================
170--- ../trunk-jpl/src/c/classes/gauss/GaussTria.h (revision 14778)
171+++ ../trunk-jpl/src/c/classes/gauss/GaussTria.h (revision 14779)
172@@ -31,6 +31,7 @@
173 GaussTria();
174 GaussTria(int order);
175 GaussTria(int index1,int index2,int order);
176+ GaussTria(int index,IssmDouble r1, IssmDouble r2,bool maintlyfloating,int order);
177 ~GaussTria();
178
179 /*Methods*/
Note: See TracBrowser for help on using the repository browser.