/[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 153 - (show annotations)
Tue Oct 25 01:51:20 2005 UTC (14 years, 1 month ago) by jgs
Original Path: trunk/esys2/finley/src/finleyC/Mesh_rec8.c
File MIME type: text/plain
File size: 10284 byte(s)
Merge of development branch dev-02 back to main trunk on 2005-10-25

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26