/[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 730 - (show annotations)
Mon May 15 04:03:49 2006 UTC (13 years, 4 months ago) by bcumming
File MIME type: text/plain
File size: 10069 byte(s)


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26