/[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 150 - (hide annotations)
Thu Sep 15 03:44:45 2005 UTC (14 years, 1 month ago) by jgs
Original Path: trunk/esys2/finley/src/finleyC/Mesh_rec4.c
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 jgs 150 /*
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 jgs 82 /**************************************************************/
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 jgs 123 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 jgs 82 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 jgs 150 if (! Finley_noError()) return NULL;
64 jgs 82
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 jgs 150 if (! Finley_noError()) {
75 jgs 82 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 jgs 150 if (! Finley_noError()) {
85 jgs 82 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 jgs 123 out->Elements->minColor=0;
136     out->Elements->maxColor=COLOR_MOD(0)+3*COLOR_MOD(0);
137 jgs 82
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 jgs 123 out->FaceElements->minColor=0;
246     out->FaceElements->maxColor=7;
247 jgs 82
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 jgs 149 #ifdef Finley_TRACE
259 jgs 82 printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);
260 jgs 149 #endif
261 jgs 82
262 jgs 150 if (! Finley_noError()) {
263 jgs 82 Finley_Mesh_dealloc(out);
264     return NULL;
265     }
266     return out;
267     }
268    
269     /*
270 jgs 149 * Merge of development branch dev-02 back to main trunk on 2005-09-01
271     *
272 jgs 150 * 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 jgs 149 * 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 jgs 123 * Revision 1.2 2005/07/08 04:07:54 jgs
279     * Merge of development branch back to main trunk on 2005-07-08
280 jgs 82 *
281 jgs 123 * 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 jgs 82 * 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