source: issm/oecreview/Archive/17984-18295/ISSM-18170-18171.diff@ 18296

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

Added 17984-18295
File size: 14.9 KB
RevLine 
[18296]1Index: ../trunk-jpl/src/c/classes/Elements/PentaRef.cpp
2===================================================================
3--- ../trunk-jpl/src/c/classes/Elements/PentaRef.cpp (revision 18170)
4+++ ../trunk-jpl/src/c/classes/Elements/PentaRef.cpp (revision 18171)
5@@ -23,6 +23,7 @@
6 #define NUMNODESP1xP3 12
7 #define NUMNODESP2xP1 12
8 #define NUMNODESP2 18
9+#define NUMNODESP2b 19
10 #define NUMNODESP2xP4 30
11
12 /*Object constructors and destructor*/
13@@ -237,6 +238,32 @@
14 basis[16]=4.*gauss->coord3*gauss->coord1*(1.-zeta*zeta);
15 basis[17]=4.*gauss->coord1*gauss->coord2*(1.-zeta*zeta);
16 return;
17+ case P2bubbleEnum: case P2bubblecondensedEnum:
18+ /*Corner nodes*/
19+ basis[ 0]=gauss->coord1*(2.*gauss->coord1-1.)*zeta*(zeta-1.)/2.;
20+ basis[ 1]=gauss->coord2*(2.*gauss->coord2-1.)*zeta*(zeta-1.)/2.;
21+ basis[ 2]=gauss->coord3*(2.*gauss->coord3-1.)*zeta*(zeta-1.)/2.;
22+ basis[ 3]=gauss->coord1*(2.*gauss->coord1-1.)*zeta*(zeta+1.)/2.;
23+ basis[ 4]=gauss->coord2*(2.*gauss->coord2-1.)*zeta*(zeta+1.)/2.;
24+ basis[ 5]=gauss->coord3*(2.*gauss->coord3-1.)*zeta*(zeta+1.)/2.;
25+ /*mid-sides of quads*/
26+ basis[ 6]=gauss->coord1*(2.*gauss->coord1-1.)*(1.-zeta*zeta);
27+ basis[ 7]=gauss->coord2*(2.*gauss->coord2-1.)*(1.-zeta*zeta);
28+ basis[ 8]=gauss->coord3*(2.*gauss->coord3-1.)*(1.-zeta*zeta);
29+ /*mid-sides of triangles*/
30+ basis[ 9]=4.*gauss->coord3*gauss->coord2*zeta*(zeta-1.)/2.;
31+ basis[10]=4.*gauss->coord3*gauss->coord1*zeta*(zeta-1.)/2.;
32+ basis[11]=4.*gauss->coord1*gauss->coord2*zeta*(zeta-1.)/2.;
33+ basis[12]=4.*gauss->coord3*gauss->coord2*zeta*(zeta+1.)/2.;
34+ basis[13]=4.*gauss->coord3*gauss->coord1*zeta*(zeta+1.)/2.;
35+ basis[14]=4.*gauss->coord1*gauss->coord2*zeta*(zeta+1.)/2.;
36+ /*quad faces*/
37+ basis[15]=4.*gauss->coord3*gauss->coord2*(1.-zeta*zeta);
38+ basis[16]=4.*gauss->coord3*gauss->coord1*(1.-zeta*zeta);
39+ basis[17]=4.*gauss->coord1*gauss->coord2*(1.-zeta*zeta);
40+ /*bubble*/
41+ basis[18]=27.*gauss->coord1*gauss->coord2*gauss->coord3*(1.+zeta)*(1.-zeta);
42+ return;
43 case P2xP4Enum :
44 /*Corner nodes*/
45 basis[ 0]=gauss->coord1*(2.*gauss->coord1-1.)*(2./3.)*(zeta-1.)*(zeta-0.5 )*(zeta)*(zeta+0.5);
46@@ -572,6 +599,88 @@
47 dbasis[NUMNODESP2*1+17] = (1.-zeta*zeta)*(-2.*SQRT3/3.*(gauss->coord1+gauss->coord2));
48 dbasis[NUMNODESP2*2+17] = -2.*zeta*4.*gauss->coord1*gauss->coord2;
49 return;
50+ case P2bubbleEnum: case P2bubblecondensedEnum:
51+ /*Nodal function 1*/
52+ dbasis[NUMNODESP2b*0+0 ] = .5*zeta*(zeta-1.)*(-2.*gauss->coord1 + 0.5);
53+ dbasis[NUMNODESP2b*1+0 ] = .5*zeta*(zeta-1.)*(-2.*SQRT3/3.*gauss->coord1 + SQRT3/6.);
54+ dbasis[NUMNODESP2b*2+0 ] = .5*(2.*zeta-1.)*gauss->coord1*(2.*gauss->coord1-1.);
55+ /*Nodal function 2*/
56+ dbasis[NUMNODESP2b*0+1 ] = .5*zeta*(zeta-1.)*(+2.*gauss->coord2 - 0.5);
57+ dbasis[NUMNODESP2b*1+1 ] = .5*zeta*(zeta-1.)*(-2.*SQRT3/3.*gauss->coord2 + SQRT3/6.);
58+ dbasis[NUMNODESP2b*2+1 ] = .5*(2.*zeta-1.)*gauss->coord2*(2.*gauss->coord2-1.);
59+ /*Nodal function 3*/
60+ dbasis[NUMNODESP2b*0+2 ] = 0.;
61+ dbasis[NUMNODESP2b*1+2 ] = .5*zeta*(zeta-1.)*(4.*SQRT3/3.*gauss->coord3 - SQRT3/3.);
62+ dbasis[NUMNODESP2b*2+2 ] = .5*(2.*zeta-1.)*gauss->coord3*(2.*gauss->coord3-1.);
63+ /*Nodal function 4*/
64+ dbasis[NUMNODESP2b*0+3 ] = .5*zeta*(zeta+1.)*(-2.*gauss->coord1 + 0.5);
65+ dbasis[NUMNODESP2b*1+3 ] = .5*zeta*(zeta+1.)*(-2.*SQRT3/3.*gauss->coord1 + SQRT3/6.);
66+ dbasis[NUMNODESP2b*2+3 ] = .5*(2.*zeta+1.)*gauss->coord1*(2.*gauss->coord1-1.);
67+ /*Nodal function 5*/
68+ dbasis[NUMNODESP2b*0+4 ] = .5*zeta*(zeta+1.)*(+2.*gauss->coord2 - 0.5);
69+ dbasis[NUMNODESP2b*1+4 ] = .5*zeta*(zeta+1.)*(-2.*SQRT3/3.*gauss->coord2 + SQRT3/6.);
70+ dbasis[NUMNODESP2b*2+4 ] = .5*(2.*zeta+1.)*gauss->coord2*(2.*gauss->coord2-1.);
71+ /*Nodal function 6*/
72+ dbasis[NUMNODESP2b*0+5 ] = 0.;
73+ dbasis[NUMNODESP2b*1+5 ] = .5*zeta*(zeta+1.)*(4.*SQRT3/3.*gauss->coord3 - SQRT3/3.);
74+ dbasis[NUMNODESP2b*2+5 ] = .5*(2.*zeta+1.)*gauss->coord3*(2.*gauss->coord3-1.);
75+
76+ /*Nodal function 7*/
77+ dbasis[NUMNODESP2b*0+6 ] = (-2.*gauss->coord1 + 0.5)*(1.-zeta*zeta);
78+ dbasis[NUMNODESP2b*1+6 ] = (-2.*SQRT3/3.*gauss->coord1 + SQRT3/6.)*(1.-zeta*zeta);
79+ dbasis[NUMNODESP2b*2+6 ] = -2.*zeta*gauss->coord1*(2.*gauss->coord1-1.);
80+ /*Nodal function 8*/
81+ dbasis[NUMNODESP2b*0+7 ] = (+2.*gauss->coord2 - 0.5)*(1.-zeta*zeta);
82+ dbasis[NUMNODESP2b*1+7 ] = (-2.*SQRT3/3.*gauss->coord2 + SQRT3/6.)*(1.-zeta*zeta);
83+ dbasis[NUMNODESP2b*2+7 ] = -2.*zeta*gauss->coord2*(2.*gauss->coord2-1.);
84+ /*Nodal function 9*/
85+ dbasis[NUMNODESP2b*0+8 ] = 0.;
86+ dbasis[NUMNODESP2b*1+8 ] = (+4.*SQRT3/3.*gauss->coord3 - SQRT3/3.)*(1.-zeta*zeta);
87+ dbasis[NUMNODESP2b*2+8 ] = -2.*zeta*gauss->coord3*(2.*gauss->coord3-1.);
88+
89+ /*Nodal function 10*/
90+ dbasis[NUMNODESP2b*0+9 ] = zeta*(zeta-1.)*gauss->coord3;
91+ dbasis[NUMNODESP2b*1+9 ] = .5*zeta*(zeta-1.)*(+4.*SQRT3/3.*gauss->coord2 - 2.*SQRT3/3.*gauss->coord3);
92+ dbasis[NUMNODESP2b*2+9 ] = 2.*gauss->coord3*gauss->coord2*(2.*zeta-1.);
93+ /*Nodal function 11*/
94+ dbasis[NUMNODESP2b*0+10] = -zeta*(zeta-1.)*gauss->coord3;
95+ dbasis[NUMNODESP2b*1+10] = .5*zeta*(zeta-1.)*(+4.*SQRT3/3.*gauss->coord1 - 2.*SQRT3/3.*gauss->coord3);
96+ dbasis[NUMNODESP2b*2+10] = 2.*gauss->coord3*gauss->coord1*(2.*zeta-1.);
97+ /*Nodal function 12*/
98+ dbasis[NUMNODESP2b*0+11] = zeta*(zeta-1.)*(gauss->coord1-gauss->coord2);
99+ dbasis[NUMNODESP2b*1+11] = .5*zeta*(zeta-1.)*(-2.*SQRT3/3.*(gauss->coord1+gauss->coord2));
100+ dbasis[NUMNODESP2b*2+11] = 2.*gauss->coord1*gauss->coord2*(2.*zeta-1.);
101+ /*Nodal function 13*/
102+ dbasis[NUMNODESP2b*0+12] = zeta*(zeta+1.)*gauss->coord3;
103+ dbasis[NUMNODESP2b*1+12] = .5*zeta*(zeta+1.)*(+4.*SQRT3/3.*gauss->coord2 - 2.*SQRT3/3.*gauss->coord3);
104+ dbasis[NUMNODESP2b*2+12] = 2.*gauss->coord3*gauss->coord2*(2.*zeta+1.);
105+ /*Nodal function 14*/
106+ dbasis[NUMNODESP2b*0+13] = -zeta*(zeta+1.)*gauss->coord3;
107+ dbasis[NUMNODESP2b*1+13] = .5*zeta*(zeta+1.)*(+4.*SQRT3/3.*gauss->coord1 - 2.*SQRT3/3.*gauss->coord3);
108+ dbasis[NUMNODESP2b*2+13] = 2.*gauss->coord3*gauss->coord1*(2.*zeta+1.);
109+ /*Nodal function 15*/
110+ dbasis[NUMNODESP2b*0+14] = zeta*(zeta+1.)*(gauss->coord1-gauss->coord2);
111+ dbasis[NUMNODESP2b*1+14] = .5*zeta*(zeta+1.)*(-2.*SQRT3/3.*(gauss->coord1+gauss->coord2));
112+ dbasis[NUMNODESP2b*2+14] = 2.*gauss->coord1*gauss->coord2*(2.*zeta+1.);
113+
114+ /*Nodal function 16*/
115+ dbasis[NUMNODESP2b*0+15] = 2.*gauss->coord3*(1.-zeta*zeta);
116+ dbasis[NUMNODESP2b*1+15] = (1.-zeta*zeta)*(+4.*SQRT3/3.*gauss->coord2 - 2.*SQRT3/3.*gauss->coord3);
117+ dbasis[NUMNODESP2b*2+15] = -2.*zeta*4.*gauss->coord3*gauss->coord2;
118+ /*Nodal function 17*/
119+ dbasis[NUMNODESP2b*0+16] = -2.*gauss->coord3*(1.-zeta*zeta);
120+ dbasis[NUMNODESP2b*1+16] = (1.-zeta*zeta)*(+4.*SQRT3/3.*gauss->coord1 - 2.*SQRT3/3.*gauss->coord3);
121+ dbasis[NUMNODESP2b*2+16] = -2.*zeta*4.*gauss->coord3*gauss->coord1;
122+ /*Nodal function 18*/
123+ dbasis[NUMNODESP2b*0+17] = 2.*(gauss->coord1-gauss->coord2)*(1.-zeta*zeta);
124+ dbasis[NUMNODESP2b*1+17] = (1.-zeta*zeta)*(-2.*SQRT3/3.*(gauss->coord1+gauss->coord2));
125+ dbasis[NUMNODESP2b*2+17] = -2.*zeta*4.*gauss->coord1*gauss->coord2;
126+
127+ /*Nodal function 19*/
128+ dbasis[NUMNODESP2b*0+18] = 27.*(1.+zeta)*(1.-zeta)*(-.5*gauss->coord2*gauss->coord3 + .5*gauss->coord1*gauss->coord3);
129+ dbasis[NUMNODESP2b*1+18] = 27.*(1.+zeta)*(1.-zeta)*SQRT3*(-1./6.*gauss->coord2*gauss->coord3 - 1./6.*gauss->coord1*gauss->coord3 +1./3.*gauss->coord1*gauss->coord2);
130+ dbasis[NUMNODESP2b*2+18] = -54*gauss->coord1*gauss->coord2*gauss->coord3*zeta;
131+ return;
132 case P2xP4Enum :
133 /*Nodal function 1*/
134 dbasis[NUMNODESP2xP4*0+0 ] = (-2* gauss->coord1 + 0.5 ) *(2./3.) *(zeta - 1.)*(zeta - 0.5)*(zeta)*(zeta+0.5) ;
135@@ -840,6 +949,8 @@
136 case P1bubbleEnum: return NUMNODESP1b;
137 case P1bubblecondensedEnum: return NUMNODESP1b;
138 case P2Enum: return NUMNODESP2;
139+ case P2bubbleEnum: return NUMNODESP2b;
140+ case P2bubblecondensedEnum: return NUMNODESP2b;
141 case P2xP1Enum: return NUMNODESP2xP1;
142 case P1xP2Enum: return NUMNODESP1xP2;
143 case P1P1Enum: return NUMNODESP1*2;
144Index: ../trunk-jpl/src/c/classes/Elements/Penta.cpp
145===================================================================
146--- ../trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 18170)
147+++ ../trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 18171)
148@@ -2233,6 +2233,29 @@
149 penta_node_ids[16]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+iomodel->elementtofaceconnectivity[5*index+3]+1;
150 penta_node_ids[17]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+iomodel->elementtofaceconnectivity[5*index+4]+1;
151 break;
152+ case P2bubbleEnum: case P2bubblecondensedEnum:
153+ numnodes = 19;
154+ penta_node_ids = xNew<int>(numnodes);
155+ penta_node_ids[ 0]=iomodel->nodecounter+iomodel->elements[6*index+0];
156+ penta_node_ids[ 1]=iomodel->nodecounter+iomodel->elements[6*index+1];
157+ penta_node_ids[ 2]=iomodel->nodecounter+iomodel->elements[6*index+2];
158+ penta_node_ids[ 3]=iomodel->nodecounter+iomodel->elements[6*index+3];
159+ penta_node_ids[ 4]=iomodel->nodecounter+iomodel->elements[6*index+4];
160+ penta_node_ids[ 5]=iomodel->nodecounter+iomodel->elements[6*index+5];
161+ penta_node_ids[ 6]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+0]+1;
162+ penta_node_ids[ 7]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+1]+1;
163+ penta_node_ids[ 8]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+2]+1;
164+ penta_node_ids[ 9]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+3]+1;
165+ penta_node_ids[10]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+4]+1;
166+ penta_node_ids[11]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+5]+1;
167+ penta_node_ids[12]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+6]+1;
168+ penta_node_ids[13]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+7]+1;
169+ penta_node_ids[14]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+8]+1;
170+ penta_node_ids[15]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+iomodel->elementtofaceconnectivity[5*index+2]+1;
171+ penta_node_ids[16]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+iomodel->elementtofaceconnectivity[5*index+3]+1;
172+ penta_node_ids[17]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+iomodel->elementtofaceconnectivity[5*index+4]+1;
173+ penta_node_ids[18]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberoffaces+index+1;
174+ break;
175 case P1P1Enum: case P1P1GLSEnum:
176 numnodes = 12;
177 penta_node_ids = xNew<int>(numnodes);
178Index: ../trunk-jpl/src/c/classes/gauss/GaussPenta.cpp
179===================================================================
180--- ../trunk-jpl/src/c/classes/gauss/GaussPenta.cpp (revision 18170)
181+++ ../trunk-jpl/src/c/classes/gauss/GaussPenta.cpp (revision 18171)
182@@ -646,6 +646,34 @@
183 default: _error_("node index should be in [0 17]");
184 }
185 break;
186+ case P2bubbleEnum: case P2bubblecondensedEnum:
187+ switch(iv){
188+ case 0: coord1=1.; coord2=0.; coord3=0.; coord4=-1.; break;
189+ case 1: coord1=0.; coord2=1.; coord3=0.; coord4=-1.; break;
190+ case 2: coord1=0.; coord2=0.; coord3=1.; coord4=-1.; break;
191+ case 3: coord1=1.; coord2=0.; coord3=0.; coord4=+1.; break;
192+ case 4: coord1=0.; coord2=1.; coord3=0.; coord4=+1.; break;
193+ case 5: coord1=0.; coord2=0.; coord3=1.; coord4=+1.; break;
194+
195+ case 6: coord1=1.; coord2=0.; coord3=0.; coord4=0.; break;
196+ case 7: coord1=0.; coord2=1.; coord3=0.; coord4=0.; break;
197+ case 8: coord1=0.; coord2=0.; coord3=1.; coord4=0.; break;
198+
199+ case 9: coord1=0.; coord2=.5; coord3=.5; coord4=-1.;break;
200+ case 10: coord1=.5; coord2=0.; coord3=.5; coord4=-1.;break;
201+ case 11: coord1=.5; coord2=.5; coord3=0.; coord4=-1.;break;
202+ case 12: coord1=0.; coord2=.5; coord3=.5; coord4=+1.;break;
203+ case 13: coord1=.5; coord2=0.; coord3=.5; coord4=+1.;break;
204+ case 14: coord1=.5; coord2=.5; coord3=0.; coord4=+1.;break;
205+
206+ case 15: coord1=0.; coord2=.5; coord3=.5; coord4=0.;break;
207+ case 16: coord1=.5; coord2=0.; coord3=.5; coord4=0.;break;
208+ case 17: coord1=.5; coord2=.5; coord3=0.; coord4=0.;break;
209+
210+ case 18: coord1=1./3.; coord2=1./3.; coord3=1./3.; coord4=0.; break;
211+ default: _error_("node index should be in [0 18]");
212+ }
213+ break;
214 case P2xP4Enum:
215 switch(iv){
216 case 0: coord1=1.; coord2=0.; coord3=0.; coord4=-1.; break;
217Index: ../trunk-jpl/src/m/classes/flowequation.py
218===================================================================
219--- ../trunk-jpl/src/m/classes/flowequation.py (revision 18170)
220+++ ../trunk-jpl/src/m/classes/flowequation.py (revision 18171)
221@@ -76,7 +76,7 @@
222 md = checkfield(md,'fieldname','flowequation.isHO','numel',[1],'values',[0,1])
223 md = checkfield(md,'fieldname','flowequation.isFS','numel',[1],'values',[0,1])
224 md = checkfield(md,'fieldname','flowequation.fe_SSA','values',['P1','P1bubble','P1bubblecondensed','P2','P2bubble'])
225- md = checkfield(md,'fieldname','flowequation.fe_HO' ,'values',['P1','P1bubble','P1bubblecondensed','P1xP2','P2xP1','P2','P1xP3','P2xP4'])
226+ md = checkfield(md,'fieldname','flowequation.fe_HO' ,'values',['P1','P1bubble','P1bubblecondensed','P1xP2','P2xP1','P2','P2bubble','P1xP3','P2xP4'])
227 md = checkfield(md,'fieldname','flowequation.fe_FS' ,'values',['P1P1','P1P1GLS','MINIcondensed','MINI','TaylorHood','XTaylorHood','OneLayerP4z'])
228 md = checkfield(md,'fieldname','flowequation.borderSSA','size',[md.mesh.numberofvertices],'values',[0,1])
229 md = checkfield(md,'fieldname','flowequation.borderHO','size',[md.mesh.numberofvertices],'values',[0,1])
230Index: ../trunk-jpl/src/m/classes/flowequation.m
231===================================================================
232--- ../trunk-jpl/src/m/classes/flowequation.m (revision 18170)
233+++ ../trunk-jpl/src/m/classes/flowequation.m (revision 18171)
234@@ -132,7 +132,7 @@
235 md = checkfield(md,'fieldname','flowequation.isHO','numel',[1],'values',[0 1]);
236 md = checkfield(md,'fieldname','flowequation.isFS','numel',[1],'values',[0 1]);
237 md = checkfield(md,'fieldname','flowequation.fe_SSA','values',{'P1','P1bubble','P1bubblecondensed','P2','P2bubble'});
238- md = checkfield(md,'fieldname','flowequation.fe_HO' ,'values',{'P1','P1bubble','P1bubblecondensed','P1xP2','P2xP1','P2','P1xP3','P2xP4'});
239+ md = checkfield(md,'fieldname','flowequation.fe_HO' ,'values',{'P1','P1bubble','P1bubblecondensed','P1xP2','P2xP1','P2','P2bubble','P1xP3','P2xP4'});
240 md = checkfield(md,'fieldname','flowequation.fe_FS' ,'values',{'P1P1','P1P1GLS','MINIcondensed','MINI','TaylorHood','XTaylorHood','OneLayerP4z'});
241 md = checkfield(md,'fieldname','flowequation.XTH_r','numel',[1],'>',0.);
242 md = checkfield(md,'fieldname','flowequation.XTH_theta','numel',[1],'>=',0.,'<',0.5);
Note: See TracBrowser for help on using the repository browser.