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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 616 - (hide annotations)
Wed Mar 22 02:46:56 2006 UTC (13 years, 7 months ago) by elspeth
File MIME type: text/plain
File size: 8509 byte(s)
Copyright added to more source files.

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 first order elements (Rec4) in the rectangle */
18     /* [0,Length[0]] x [0,Length[1]]. order is the desired accuracy of the integration scheme. */
19    
20    
21     /**************************************************************/
22    
23     /* Author: gross@access.edu.au */
24     /* Version: $Id$ */
25    
26     /**************************************************************/
27    
28     #include "RectangularMesh.h"
29    
30     /**************************************************************/
31    
32 jgs 123 Finley_Mesh* Finley_RectangularMesh_Rec4(dim_t* numElements,double* Length,bool_t* periodic, index_t order,bool_t useElementsOnFace) {
33 jgs 153 dim_t N0,N1,NE0,NE1,i0,i1,totalNECount,faceNECount,NDOF0,NDOF1,NFaceElements,NUMNODES,M0,M1;
34 jgs 123 index_t k,node0;
35 jgs 82 Finley_Mesh* out;
36     char name[50];
37     double time0=Finley_timer();
38     NE0=MAX(1,numElements[0]);
39     NE1=MAX(1,numElements[1]);
40     N0=NE0+1;
41     N1=NE1+1;
42    
43 jgs 153 if (N0<=N1) {
44     M0=1;
45     M1=N0;
46     } else {
47     M0=N1;
48     M1=1;
49     }
50    
51 jgs 82 NFaceElements=0;
52     if (!periodic[0]) {
53     NDOF0=N0;
54     NFaceElements+=2*NE1;
55     } else {
56     NDOF0=N0-1;
57     }
58     if (!periodic[1]) {
59     NDOF1=N1;
60     NFaceElements+=2*NE0;
61     } else {
62     NDOF1=N1-1;
63     }
64    
65     /* allocate mesh: */
66    
67     sprintf(name,"Rectangular %d x %d mesh",N0,N1);
68     out=Finley_Mesh_alloc(name,2,order);
69 jgs 150 if (! Finley_noError()) return NULL;
70 jgs 82
71     out->Elements=Finley_ElementFile_alloc(Rec4,out->order);
72     if (useElementsOnFace) {
73     out->FaceElements=Finley_ElementFile_alloc(Rec4Face,out->order);
74     out->ContactElements=Finley_ElementFile_alloc(Rec4Face_Contact,out->order);
75     } else {
76     out->FaceElements=Finley_ElementFile_alloc(Line2,out->order);
77     out->ContactElements=Finley_ElementFile_alloc(Line2_Contact,out->order);
78     }
79     out->Points=Finley_ElementFile_alloc(Point1,out->order);
80 jgs 150 if (! Finley_noError()) {
81 jgs 82 Finley_Mesh_dealloc(out);
82     return NULL;
83     }
84    
85     /* allocate tables: */
86    
87     Finley_NodeFile_allocTable(out->Nodes,N0*N1);
88     Finley_ElementFile_allocTable(out->Elements,NE0*NE1);
89     Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);
90 jgs 150 if (! Finley_noError()) {
91 jgs 82 Finley_Mesh_dealloc(out);
92     return NULL;
93     }
94     /* set nodes: */
95 jgs 153
96     #pragma omp parallel for private(i0,i1,k)
97 jgs 82 for (i1=0;i1<N1;i1++) {
98     for (i0=0;i0<N0;i0++) {
99 jgs 153 k=M0*i0+M1*i1;
100 jgs 82 out->Nodes->Coordinates[INDEX2(0,k,2)]=DBLE(i0)/DBLE(N0-1)*Length[0];
101     out->Nodes->Coordinates[INDEX2(1,k,2)]=DBLE(i1)/DBLE(N1-1)*Length[1];
102 jgs 153 out->Nodes->Id[k]=i0+N0*i1;
103 jgs 82 out->Nodes->Tag[k]=0;
104 jgs 153 out->Nodes->degreeOfFreedom[k]=M0*(i0%NDOF0) +M1*(i1%NDOF1);
105 jgs 82 }
106     }
107     /* tags for the faces: */
108     if (!periodic[1]) {
109     for (i0=0;i0<N0;i0++) {
110 jgs 153 out->Nodes->Tag[M0*i0+M1*0]+=10;
111     out->Nodes->Tag[M0*i0+M1*(N1-1)]+=20;
112 jgs 82 }
113     }
114     if (!periodic[0]) {
115     for (i1=0;i1<N1;i1++) {
116 jgs 153 out->Nodes->Tag[M0*0+M1*i1]+=1;
117     out->Nodes->Tag[M0*(N0-1)+M1*i1]+=2;
118 jgs 82 }
119     }
120 jgs 153
121 jgs 82 /* set the elements: */
122    
123     #pragma omp parallel for private(i0,i1,k,node0)
124     for (i1=0;i1<NE1;i1++) {
125     for (i0=0;i0<NE0;i0++) {
126     k=i0+NE0*i1;
127     node0=i0+i1*N0;
128    
129     out->Elements->Id[k]=k;
130     out->Elements->Tag[k]=0;
131     out->Elements->Color[k]=COLOR_MOD(i0)+3*COLOR_MOD(i1);
132    
133     out->Elements->Nodes[INDEX2(0,k,4)]=node0;
134     out->Elements->Nodes[INDEX2(1,k,4)]=node0+1;
135     out->Elements->Nodes[INDEX2(2,k,4)]=node0+N0+1;
136     out->Elements->Nodes[INDEX2(3,k,4)]=node0+N0;
137    
138     }
139     }
140 jgs 123 out->Elements->minColor=0;
141     out->Elements->maxColor=COLOR_MOD(0)+3*COLOR_MOD(0);
142 jgs 82
143     /* face elements: */
144    
145     if (useElementsOnFace) {
146     NUMNODES=4;
147     } else {
148     NUMNODES=2;
149     }
150     totalNECount=NE0*NE1;
151     faceNECount=0;
152    
153     /* elements on boundary 010 (x2=0): */
154     if (!periodic[1]) {
155     #pragma omp parallel for private(i0,k,node0)
156     for (i0=0;i0<NE0;i0++) {
157     k=i0+faceNECount;
158     node0=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+1;
167     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0+1;
168     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0;
169     } else {
170     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
171     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+1;
172     }
173     }
174     totalNECount+=NE0;
175     faceNECount+=NE0;
176    
177     /* ** elements on boundary 020 (x2=1): */
178    
179     #pragma omp parallel for private(i0,k,node0)
180     for (i0=0;i0<NE0;i0++) {
181     k=i0+faceNECount;
182     node0=i0+(NE1-1)*N0;
183    
184     out->FaceElements->Id[k]=i0+totalNECount;
185     out->FaceElements->Tag[k]=20;
186     out->FaceElements->Color[k]=i0%2+2;
187    
188     if (useElementsOnFace) {
189     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+N0+1;
190     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0;
191     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0;
192     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+1;
193     } else {
194     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+N0+1;
195     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0;
196     }
197     }
198     totalNECount+=NE0;
199     faceNECount+=NE0;
200     }
201     if (!periodic[0]) {
202     /* ** elements on boundary 001 (x1=0): */
203     #pragma omp parallel for private(i1,k,node0)
204     for (i1=0;i1<NE1;i1++) {
205     k=i1+faceNECount;
206     node0=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+N0;
214     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0;
215     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+1;
216     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0+1;
217     } else {
218     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+N0;
219     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0;
220    
221     }
222     }
223     totalNECount+=NE1;
224     faceNECount+=NE1;
225    
226     /* ** elements on boundary 002 (x1=1): */
227    
228     #pragma omp parallel for private(i1,k,node0)
229     for (i1=0;i1<NE1;i1++) {
230     k=i1+faceNECount;
231     node0=(NE0-1)+i1*N0;
232    
233     out->FaceElements->Id[k]=i1+totalNECount;
234     out->FaceElements->Tag[k]=2;
235     out->FaceElements->Color[k]=i1%2+6;
236    
237     if (useElementsOnFace) {
238     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+1;
239     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0+1;
240     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0;
241     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0;
242     } else {
243     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+1;
244     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0+1;
245     }
246     }
247     totalNECount+=NE1;
248     faceNECount+=NE1;
249     }
250 jgs 123 out->FaceElements->minColor=0;
251     out->FaceElements->maxColor=7;
252 jgs 82
253     /* face elements done: */
254    
255     /* condense the nodes: */
256    
257     Finley_Mesh_resolveNodeIds(out);
258    
259     /* prepare mesh for further calculatuions:*/
260    
261     Finley_Mesh_prepare(out) ;
262    
263 jgs 149 #ifdef Finley_TRACE
264 jgs 82 printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);
265 jgs 149 #endif
266 jgs 82
267 jgs 150 if (! Finley_noError()) {
268 jgs 82 Finley_Mesh_dealloc(out);
269     return NULL;
270     }
271     return out;
272     }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26