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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 730 - (hide annotations)
Mon May 15 04:03:49 2006 UTC (13 years, 5 months ago) by bcumming
File MIME type: text/plain
File size: 10069 byte(s)


1 jgs 150 /*
2 elspeth 616 ************************************************************
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 jgs 150 */
12    
13 jgs 82 /**************************************************************/
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 jgs 153 dim_t N0,N1,NE0,NE1,i0,i1,totalNECount,faceNECount,NDOF0,NDOF1,NFaceElements,NUMNODES,M0,M1;
33 jgs 123 index_t k,node0;
34 bcumming 730 Finley_Mesh* out=NULL;
35 jgs 82 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 jgs 153 if (N0<=N1) {
43     M0=1;
44     M1=N0;
45     } else {
46     M0=N1;
47     M1=1;
48     }
49    
50 jgs 82 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 bcumming 730 /* TEMPFIX */
68     #ifndef PASO_MPI
69 jgs 82 out=Finley_Mesh_alloc(name,2,order);
70 bcumming 730
71 jgs 150 if (! Finley_noError()) return NULL;
72 jgs 82
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 bcumming 730 #else
83     /* TODO */
84     #endif
85 jgs 150 if (! Finley_noError()) {
86 jgs 82 Finley_Mesh_dealloc(out);
87     return NULL;
88     }
89    
90     /* allocate tables: */
91 bcumming 730 #ifndef PASO_MPI
92 jgs 82 Finley_NodeFile_allocTable(out->Nodes,N0*N1);
93 bcumming 730 #else
94     /* TODO */
95     #endif
96 jgs 82 Finley_ElementFile_allocTable(out->Elements,NE0*NE1);
97     Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);
98 jgs 150 if (! Finley_noError()) {
99 jgs 82 Finley_Mesh_dealloc(out);
100     return NULL;
101     }
102     /* set nodes: */
103 jgs 153
104     #pragma omp parallel for private(i0,i1,k)
105 jgs 82 for (i1=0;i1<N1;i1++) {
106     for (i0=0;i0<N0;i0++) {
107 jgs 153 k=M0*i0+M1*i1;
108 jgs 82 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 jgs 153 out->Nodes->Id[k]=i0+N0*i1;
111 jgs 82 out->Nodes->Tag[k]=0;
112 jgs 153 out->Nodes->degreeOfFreedom[k]=M0*(i0%NDOF0) +M1*(i1%NDOF1);
113 jgs 82 }
114     }
115     /* tags for the faces: */
116     if (!periodic[1]) {
117     for (i0=0;i0<N0;i0++) {
118 jgs 153 out->Nodes->Tag[M0*i0+M1*0]+=10;
119     out->Nodes->Tag[M0*i0+M1*(N1-1)]+=20;
120 jgs 82 }
121     }
122     if (!periodic[0]) {
123     for (i1=0;i1<N1;i1++) {
124 jgs 153 out->Nodes->Tag[M0*0+M1*i1]+=1;
125     out->Nodes->Tag[M0*(N0-1)+M1*i1]+=2;
126 jgs 82 }
127     }
128 jgs 153
129 jgs 82 /* 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 jgs 123 out->Elements->minColor=0;
153     out->Elements->maxColor=COLOR_MOD(0)+3*COLOR_MOD(0);
154 jgs 82
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 jgs 123 out->FaceElements->minColor=0;
286     out->FaceElements->maxColor=7;
287 jgs 82
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 jgs 149 #ifdef Finley_TRACE
299 jgs 82 printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);
300 jgs 149 #endif
301 jgs 82
302 jgs 150 if (! Finley_noError()) {
303 jgs 82 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