/[escript]/trunk/finley/src/Mesh_rec8.c
ViewVC logotype

Contents of /trunk/finley/src/Mesh_rec8.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 82 - (show annotations)
Tue Oct 26 06:53:54 2004 UTC (15 years ago) by jgs
Original Path: trunk/esys2/finley/src/finleyC/Mesh_rec8.c
File MIME type: text/plain
File size: 9305 byte(s)
Initial revision

1 /**************************************************************/
2
3 /* Finley: generates rectangular meshes */
4
5 /* Generates a numElements[0] x numElements[1] mesh with second order elements (Rec8) in the rectangle */
6 /* [0,Length[0]] x [0,Length[1]]. order is the desired accuracy of the integration scheme. */
7
8 /**************************************************************/
9
10 /* Copyrights by ACcESS Australia 2003 */
11 /* Author: gross@access.edu.au */
12 /* Version: $Id$ */
13
14 /**************************************************************/
15
16 #include "Common.h"
17 #include "Finley.h"
18 #include "Mesh.h"
19 #include "RectangularMesh.h"
20
21 /**************************************************************/
22
23 Finley_Mesh* Finley_RectangularMesh_Rec8(int* numElements,double* Length,int* periodic,int order,int useElementsOnFace) {
24 int N0,N1,NE0,NE1,i0,i1,k,node0,totalNECount,faceNECount,NDOF0,NDOF1,NFaceElements,NUMNODES;
25 Finley_Mesh* out;
26 char name[50];
27 double time0=Finley_timer();
28 NE0=MAX(1,numElements[0]);
29 NE1=MAX(1,numElements[1]);
30 N0=2*NE0+1;
31 N1=2*NE1+1;
32
33 NFaceElements=0;
34 if (!periodic[0]) {
35 NDOF0=N0;
36 NFaceElements+=2*NE1;
37 } else {
38 NDOF0=N0-1;
39 }
40 if (!periodic[1]) {
41 NDOF1=N1;
42 NFaceElements+=2*NE0;
43 } else {
44 NDOF1=N1-1;
45 }
46
47 /* allocate mesh: */
48
49 sprintf(name,"Rectangular %d x %d mesh",N0,N1);
50 out=Finley_Mesh_alloc(name,2,order);
51 if (Finley_ErrorCode!=NO_ERROR) return NULL;
52
53 out->Elements=Finley_ElementFile_alloc(Rec8,out->order);
54 if (useElementsOnFace) {
55 out->FaceElements=Finley_ElementFile_alloc(Rec8Face,out->order);
56 out->ContactElements=Finley_ElementFile_alloc(Rec8Face_Contact,out->order);
57 } else {
58 out->FaceElements=Finley_ElementFile_alloc(Line3,out->order);
59 out->ContactElements=Finley_ElementFile_alloc(Line3_Contact,out->order);
60 }
61 out->Points=Finley_ElementFile_alloc(Point1,out->order);
62 if (Finley_ErrorCode!=NO_ERROR) {
63 Finley_Mesh_dealloc(out);
64 return NULL;
65 }
66
67 /* allocate tables: */
68
69 Finley_NodeFile_allocTable(out->Nodes,N0*N1);
70 Finley_ElementFile_allocTable(out->Elements,NE0*NE1);
71 Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);
72 if (Finley_ErrorCode!=NO_ERROR) {
73 Finley_Mesh_dealloc(out);
74 return NULL;
75 }
76
77 /* set nodes: */
78
79 #pragma omp parallel for private(i0,i1,k)
80 for (i1=0;i1<N1;i1++) {
81 for (i0=0;i0<N0;i0++) {
82 k=i0+N0*i1;
83 out->Nodes->Coordinates[INDEX2(0,k,2)]=DBLE(i0)/DBLE(N0-1)*Length[0];
84 out->Nodes->Coordinates[INDEX2(1,k,2)]=DBLE(i1)/DBLE(N1-1)*Length[1];
85 out->Nodes->Id[k]=k;
86 out->Nodes->Tag[k]=0;
87 out->Nodes->degreeOfFreedom[k]=(i0%NDOF0) +N0*(i1%NDOF1);
88 }
89 }
90 /* tags for the faces: */
91 if (!periodic[1]) {
92 for (i0=0;i0<N0;i0++) {
93 out->Nodes->Tag[i0+N0*0]+=10;
94 out->Nodes->Tag[i0+N0*(N1-1)]+=20;
95 }
96 }
97 if (!periodic[0]) {
98 for (i1=0;i1<N1;i1++) {
99 out->Nodes->Tag[0+N0*i1]+=1;
100 out->Nodes->Tag[(N0-1)+N0*i1]+=2;
101 }
102 }
103
104 /* set the elements: */
105
106 #pragma omp parallel for private(i0,i1,k,node0)
107 for (i1=0;i1<NE1;i1++) {
108 for (i0=0;i0<NE0;i0++) {
109 k=i0+NE0*i1;
110 node0=2*i0+2*i1*N0;
111
112 out->Elements->Id[k]=k;
113 out->Elements->Tag[k]=0;
114 out->Elements->Color[k]=COLOR_MOD(i0)+3*COLOR_MOD(i1);
115
116 out->Elements->Nodes[INDEX2(0,k,8)]=node0;
117 out->Elements->Nodes[INDEX2(1,k,8)]=node0+2;
118 out->Elements->Nodes[INDEX2(2,k,8)]=node0+2*N0+2;
119 out->Elements->Nodes[INDEX2(3,k,8)]=node0+2*N0;
120 out->Elements->Nodes[INDEX2(4,k,8)]=node0+1;
121 out->Elements->Nodes[INDEX2(5,k,8)]=node0+N0+2;
122 out->Elements->Nodes[INDEX2(6,k,8)]=node0+2*N0+1;
123 out->Elements->Nodes[INDEX2(7,k,8)]=node0+N0;
124
125 }
126 }
127 out->Elements->numColors=COLOR_MOD(0)+3*COLOR_MOD(0)+1;
128
129 /* face elements: */
130 if (useElementsOnFace) {
131 NUMNODES=8;
132 } else {
133 NUMNODES=3;
134 }
135
136
137 totalNECount=NE0*NE1;
138 faceNECount=0;
139
140 if (!periodic[1]) {
141 /* ** elements on boundary 010 (x2=0): */
142
143 #pragma omp parallel for private(i0,k,node0)
144 for (i0=0;i0<NE0;i0++) {
145 k=i0+faceNECount;
146 node0=2*i0;
147
148 out->FaceElements->Id[k]=i0+totalNECount;
149 out->FaceElements->Tag[k]=10;
150 out->FaceElements->Color[k]=i0%2;
151
152 if (useElementsOnFace) {
153 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
154 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2;
155 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+2*N0+2;
156 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+2*N0;
157 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+1;
158 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N0+2;
159 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+2*N0+1;
160 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N0;
161 } else {
162 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
163 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2;
164 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+1;
165 }
166 }
167 totalNECount+=NE0;
168 faceNECount+=NE0;
169
170 /* ** elements on boundary 020 (x2=1): */
171
172 #pragma omp parallel for private(i0,k,node0)
173 for (i0=0;i0<NE0;i0++) {
174 k=i0+faceNECount;
175 node0=2*i0+2*(NE1-1)*N0;
176
177 out->FaceElements->Id[k]=i0+totalNECount;
178 out->FaceElements->Tag[k]=20;
179 out->FaceElements->Color[k]=i0%2+2;
180 if (useElementsOnFace) {
181 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+2*N0+2;
182 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2*N0;
183 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0;
184 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+2;
185 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+2*N0+1;
186 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N0;
187 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+1;
188 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N0+2;
189 } else {
190 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+2*N0+2;
191 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2*N0;
192 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+2*N0+1;
193 }
194 }
195 totalNECount+=NE0;
196 faceNECount+=NE0;
197 }
198 if (!periodic[0]) {
199 /* ** elements on boundary 001 (x1=0): */
200
201 #pragma omp parallel for private(i1,k,node0)
202 for (i1=0;i1<NE1;i1++) {
203 k=i1+faceNECount;
204 node0=2*i1*N0;
205
206 out->FaceElements->Id[k]=i1+totalNECount;
207 out->FaceElements->Tag[k]=1;
208 out->FaceElements->Color[k]=(i1%2)+4;
209
210 if (useElementsOnFace) {
211 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+2*N0;
212 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0;
213 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+2;
214 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+2*N0+2;
215 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+N0;
216 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+1;
217 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0+2;
218 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+2*N0+1;
219 } else {
220 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+2*N0;
221 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0;
222 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0;
223 }
224
225 }
226 totalNECount+=NE1;
227 faceNECount+=NE1;
228
229 /* ** elements on boundary 002 (x1=1): */
230
231 #pragma omp parallel for private(i1,k,node0)
232 for (i1=0;i1<NE1;i1++) {
233 k=i1+faceNECount;
234 node0=2*(NE0-1)+2*i1*N0;
235
236 out->FaceElements->Id[k]=i1+totalNECount;
237 out->FaceElements->Tag[k]=2;
238 out->FaceElements->Color[k]=(i1%2)+6;
239
240 if (useElementsOnFace) {
241 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+2;
242 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2*N0+2;
243 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+2*N0;
244 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0;
245 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+N0+2;
246 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+2*N0+1;
247 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0;
248 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+1;
249 } else {
250 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+2;
251 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2*N0+2;
252 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0+2;
253 }
254 }
255 totalNECount+=NE1;
256 faceNECount+=NE1;
257
258 }
259 out->FaceElements->numColors=8;
260
261 /* face elements done: */
262
263 /* condense the nodes: */
264
265 Finley_Mesh_resolveNodeIds(out);
266
267 /* prepare mesh for further calculatuions:*/
268
269 Finley_Mesh_prepare(out) ;
270
271 printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);
272
273 if (Finley_ErrorCode!=NO_ERROR) {
274 Finley_Mesh_dealloc(out);
275 return NULL;
276 }
277 return out;
278 }
279
280 /*
281 * $Log$
282 * Revision 1.1 2004/10/26 06:53:57 jgs
283 * Initial revision
284 *
285 * Revision 1.1.1.1 2004/06/24 04:00:40 johng
286 * Initial version of eys using boost-python.
287 *
288 *
289 */
290

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision

  ViewVC Help
Powered by ViewVC 1.1.26