/[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 123 - (hide annotations)
Fri Jul 8 04:08:13 2005 UTC (14 years, 3 months ago) by jgs
Original Path: trunk/esys2/finley/src/finleyC/Mesh_rec8.c
File MIME type: text/plain
File size: 9599 byte(s)
Merge of development branch back to main trunk on 2005-07-08

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26