/[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 150 - (hide annotations)
Thu Sep 15 03:44:45 2005 UTC (13 years, 10 months ago) by jgs
Original Path: trunk/esys2/finley/src/finleyC/Mesh_rec8.c
File MIME type: text/plain
File size: 10754 byte(s)
Merge of development branch dev-02 back to main trunk on 2005-09-15

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26