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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 82 - (hide annotations)
Tue Oct 26 06:53:54 2004 UTC (15 years, 1 month ago) by jgs
File MIME type: text/plain
File size: 9305 byte(s)
Initial revision

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26