/[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 153 - (hide 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 jgs 150 /*
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 jgs 82 /**************************************************************/
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 jgs 153 dim_t N0,N1,NE0,NE1,i0,i1,totalNECount,faceNECount,NDOF0,NDOF1,NFaceElements,NUMNODES,M0,M1;
35 jgs 123 index_t k,node0;
36 jgs 82 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 jgs 153 if (N0<=N1) {
45     M0=1;
46     M1=N0;
47     } else {
48     M0=N1;
49     M1=1;
50     }
51    
52 jgs 82 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 jgs 150 if (! Finley_noError()) return NULL;
71 jgs 82
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 jgs 150 if (! Finley_noError()) {
82 jgs 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 jgs 150 if (! Finley_noError()) {
92 jgs 82 Finley_Mesh_dealloc(out);
93     return NULL;
94     }
95     /* set nodes: */
96 jgs 153
97     #pragma omp parallel for private(i0,i1,k)
98 jgs 82 for (i1=0;i1<N1;i1++) {
99     for (i0=0;i0<N0;i0++) {
100 jgs 153 k=M0*i0+M1*i1;
101 jgs 82 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 jgs 153 out->Nodes->Id[k]=i0+N0*i1;
104 jgs 82 out->Nodes->Tag[k]=0;
105 jgs 153 out->Nodes->degreeOfFreedom[k]=M0*(i0%NDOF0) +M1*(i1%NDOF1);
106 jgs 82 }
107     }
108     /* tags for the faces: */
109     if (!periodic[1]) {
110     for (i0=0;i0<N0;i0++) {
111 jgs 153 out->Nodes->Tag[M0*i0+M1*0]+=10;
112     out->Nodes->Tag[M0*i0+M1*(N1-1)]+=20;
113 jgs 82 }
114     }
115     if (!periodic[0]) {
116     for (i1=0;i1<N1;i1++) {
117 jgs 153 out->Nodes->Tag[M0*0+M1*i1]+=1;
118     out->Nodes->Tag[M0*(N0-1)+M1*i1]+=2;
119 jgs 82 }
120     }
121 jgs 153
122 jgs 82 /* 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 jgs 123 out->Elements->minColor=0;
146     out->Elements->maxColor=COLOR_MOD(0)+3*COLOR_MOD(0);
147 jgs 82
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 jgs 123 out->FaceElements->minColor=0;
279     out->FaceElements->maxColor=7;
280 jgs 82
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 jgs 149 #ifdef Finley_TRACE
292 jgs 82 printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);
293 jgs 149 #endif
294 jgs 82
295 jgs 150 if (! Finley_noError()) {
296 jgs 82 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