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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 123 - (show annotations)
Fri Jul 8 04:08:13 2005 UTC (14 years, 5 months ago) by jgs
File MIME type: text/plain
File size: 9599 byte(s)
Merge of development branch back to main trunk on 2005-07-08

1 /**************************************************************/
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 dim_t N0,N1,NE0,NE1,i0,i1,totalNECount,faceNECount,NDOF0,NDOF1,NFaceElements,NUMNODES;
25 index_t k,node0;
26 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 out->Elements->minColor=0;
129 out->Elements->maxColor=COLOR_MOD(0)+3*COLOR_MOD(0);
130
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 out->FaceElements->minColor=0;
262 out->FaceElements->maxColor=7;
263
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 * Revision 1.2 2005/07/08 04:07:54 jgs
286 * Merge of development branch back to main trunk on 2005-07-08
287 *
288 * 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 * 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