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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 150 - (show annotations)
Thu Sep 15 03:44:45 2005 UTC (14 years, 2 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 /*
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 /**************************************************************/
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 dim_t N0,N1,NE0,NE1,i0,i1,totalNECount,faceNECount,NDOF0,NDOF1,NFaceElements,NUMNODES;
35 index_t k,node0;
36 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 if (! Finley_noError()) return NULL;
63
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 if (! Finley_noError()) {
74 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 if (! Finley_noError()) {
84 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 out->Elements->minColor=0;
139 out->Elements->maxColor=COLOR_MOD(0)+3*COLOR_MOD(0);
140
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 out->FaceElements->minColor=0;
272 out->FaceElements->maxColor=7;
273
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 #ifdef Finley_TRACE
285 printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);
286 #endif
287
288 if (! Finley_noError()) {
289 Finley_Mesh_dealloc(out);
290 return NULL;
291 }
292 return out;
293 }
294
295 /*
296 * 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 * 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 * 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 * Revision 1.2 2005/07/08 04:07:54 jgs
306 * Merge of development branch back to main trunk on 2005-07-08
307 *
308 * 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 * 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