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

Contents of /trunk/esys2/finley/src/finleyC/Mesh_rec4.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
File MIME type: text/plain
File size: 9264 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 first order elements (Rec4) in the rectangle */
20 /* [0,Length[0]] x [0,Length[1]]. order is the desired accuracy of the integration scheme. */
21
22
23 /**************************************************************/
24
25 /* Author: gross@access.edu.au */
26 /* Version: $Id$ */
27
28 /**************************************************************/
29
30 #include "RectangularMesh.h"
31
32 /**************************************************************/
33
34 Finley_Mesh* Finley_RectangularMesh_Rec4(dim_t* numElements,double* Length,bool_t* periodic, index_t order,bool_t useElementsOnFace) {
35 dim_t N0,N1,NE0,NE1,i0,i1,totalNECount,faceNECount,NDOF0,NDOF1,NFaceElements,NUMNODES;
36 index_t k,node0;
37 Finley_Mesh* out;
38 char name[50];
39 double time0=Finley_timer();
40 NE0=MAX(1,numElements[0]);
41 NE1=MAX(1,numElements[1]);
42 N0=NE0+1;
43 N1=NE1+1;
44
45 NFaceElements=0;
46 if (!periodic[0]) {
47 NDOF0=N0;
48 NFaceElements+=2*NE1;
49 } else {
50 NDOF0=N0-1;
51 }
52 if (!periodic[1]) {
53 NDOF1=N1;
54 NFaceElements+=2*NE0;
55 } else {
56 NDOF1=N1-1;
57 }
58
59 /* allocate mesh: */
60
61 sprintf(name,"Rectangular %d x %d mesh",N0,N1);
62 out=Finley_Mesh_alloc(name,2,order);
63 if (! Finley_noError()) return NULL;
64
65 out->Elements=Finley_ElementFile_alloc(Rec4,out->order);
66 if (useElementsOnFace) {
67 out->FaceElements=Finley_ElementFile_alloc(Rec4Face,out->order);
68 out->ContactElements=Finley_ElementFile_alloc(Rec4Face_Contact,out->order);
69 } else {
70 out->FaceElements=Finley_ElementFile_alloc(Line2,out->order);
71 out->ContactElements=Finley_ElementFile_alloc(Line2_Contact,out->order);
72 }
73 out->Points=Finley_ElementFile_alloc(Point1,out->order);
74 if (! Finley_noError()) {
75 Finley_Mesh_dealloc(out);
76 return NULL;
77 }
78
79 /* allocate tables: */
80
81 Finley_NodeFile_allocTable(out->Nodes,N0*N1);
82 Finley_ElementFile_allocTable(out->Elements,NE0*NE1);
83 Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);
84 if (! Finley_noError()) {
85 Finley_Mesh_dealloc(out);
86 return NULL;
87 }
88
89 /* set nodes: */
90
91 #pragma omp parallel for private(i0,i1,k)
92 for (i1=0;i1<N1;i1++) {
93 for (i0=0;i0<N0;i0++) {
94 k=i0+N0*i1;
95 out->Nodes->Coordinates[INDEX2(0,k,2)]=DBLE(i0)/DBLE(N0-1)*Length[0];
96 out->Nodes->Coordinates[INDEX2(1,k,2)]=DBLE(i1)/DBLE(N1-1)*Length[1];
97 out->Nodes->Id[k]=k;
98 out->Nodes->Tag[k]=0;
99 out->Nodes->degreeOfFreedom[k]=(i0%NDOF0) +N0*(i1%NDOF1);
100 }
101 }
102 /* tags for the faces: */
103 if (!periodic[1]) {
104 for (i0=0;i0<N0;i0++) {
105 out->Nodes->Tag[i0+N0*0]+=10;
106 out->Nodes->Tag[i0+N0*(N1-1)]+=20;
107 }
108 }
109 if (!periodic[0]) {
110 for (i1=0;i1<N1;i1++) {
111 out->Nodes->Tag[0+N0*i1]+=1;
112 out->Nodes->Tag[(N0-1)+N0*i1]+=2;
113 }
114 }
115
116 /* set the elements: */
117
118 #pragma omp parallel for private(i0,i1,k,node0)
119 for (i1=0;i1<NE1;i1++) {
120 for (i0=0;i0<NE0;i0++) {
121 k=i0+NE0*i1;
122 node0=i0+i1*N0;
123
124 out->Elements->Id[k]=k;
125 out->Elements->Tag[k]=0;
126 out->Elements->Color[k]=COLOR_MOD(i0)+3*COLOR_MOD(i1);
127
128 out->Elements->Nodes[INDEX2(0,k,4)]=node0;
129 out->Elements->Nodes[INDEX2(1,k,4)]=node0+1;
130 out->Elements->Nodes[INDEX2(2,k,4)]=node0+N0+1;
131 out->Elements->Nodes[INDEX2(3,k,4)]=node0+N0;
132
133 }
134 }
135 out->Elements->minColor=0;
136 out->Elements->maxColor=COLOR_MOD(0)+3*COLOR_MOD(0);
137
138 /* face elements: */
139
140 if (useElementsOnFace) {
141 NUMNODES=4;
142 } else {
143 NUMNODES=2;
144 }
145 totalNECount=NE0*NE1;
146 faceNECount=0;
147
148 /* elements on boundary 010 (x2=0): */
149 if (!periodic[1]) {
150 #pragma omp parallel for private(i0,k,node0)
151 for (i0=0;i0<NE0;i0++) {
152 k=i0+faceNECount;
153 node0=i0;
154
155 out->FaceElements->Id[k]=i0+totalNECount;
156 out->FaceElements->Tag[k]=10;
157 out->FaceElements->Color[k]=i0%2;
158
159 if (useElementsOnFace) {
160 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
161 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+1;
162 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0+1;
163 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0;
164 } else {
165 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
166 out->FaceElements->Nodes[INDEX2(1,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=i0+(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
183 if (useElementsOnFace) {
184 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+N0+1;
185 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0;
186 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0;
187 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+1;
188 } else {
189 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+N0+1;
190 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0;
191 }
192 }
193 totalNECount+=NE0;
194 faceNECount+=NE0;
195 }
196 if (!periodic[0]) {
197 /* ** elements on boundary 001 (x1=0): */
198 #pragma omp parallel for private(i1,k,node0)
199 for (i1=0;i1<NE1;i1++) {
200 k=i1+faceNECount;
201 node0=i1*N0;
202
203 out->FaceElements->Id[k]=i1+totalNECount;
204 out->FaceElements->Tag[k]=1;
205 out->FaceElements->Color[k]=i1%2+4;
206
207 if (useElementsOnFace) {
208 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+N0;
209 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0;
210 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+1;
211 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0+1;
212 } else {
213 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+N0;
214 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0;
215
216 }
217 }
218 totalNECount+=NE1;
219 faceNECount+=NE1;
220
221 /* ** elements on boundary 002 (x1=1): */
222
223 #pragma omp parallel for private(i1,k,node0)
224 for (i1=0;i1<NE1;i1++) {
225 k=i1+faceNECount;
226 node0=(NE0-1)+i1*N0;
227
228 out->FaceElements->Id[k]=i1+totalNECount;
229 out->FaceElements->Tag[k]=2;
230 out->FaceElements->Color[k]=i1%2+6;
231
232 if (useElementsOnFace) {
233 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+1;
234 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0+1;
235 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0;
236 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0;
237 } else {
238 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+1;
239 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0+1;
240 }
241 }
242 totalNECount+=NE1;
243 faceNECount+=NE1;
244 }
245 out->FaceElements->minColor=0;
246 out->FaceElements->maxColor=7;
247
248 /* face elements done: */
249
250 /* condense the nodes: */
251
252 Finley_Mesh_resolveNodeIds(out);
253
254 /* prepare mesh for further calculatuions:*/
255
256 Finley_Mesh_prepare(out) ;
257
258 #ifdef Finley_TRACE
259 printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);
260 #endif
261
262 if (! Finley_noError()) {
263 Finley_Mesh_dealloc(out);
264 return NULL;
265 }
266 return out;
267 }
268
269 /*
270 * Merge of development branch dev-02 back to main trunk on 2005-09-01
271 *
272 * Revision 1.2.2.2 2005/09/07 06:26:19 gross
273 * the solver from finley are put into the standalone package paso now
274 *
275 * Revision 1.2.2.1 2005/08/24 02:02:18 gross
276 * timing output switched off. solver output can be swiched through getSolution(verbose=True) now.
277 *
278 * Revision 1.2 2005/07/08 04:07:54 jgs
279 * Merge of development branch back to main trunk on 2005-07-08
280 *
281 * Revision 1.1.1.1.2.1 2005/06/29 02:34:53 gross
282 * some changes towards 64 integers in finley
283 *
284 * Revision 1.1.1.1 2004/10/26 06:53:57 jgs
285 * initial import of project esys2
286 *
287 * Revision 1.1.1.1 2004/06/24 04:00:40 johng
288 * Initial version of eys using boost-python.
289 *
290 *
291 */
292

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26