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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 153 - (hide annotations)
Tue Oct 25 01:51:20 2005 UTC (14 years ago) by jgs
Original Path: trunk/esys2/finley/src/finleyC/Mesh_rec4.c
File MIME type: text/plain
File size: 8831 byte(s)
Merge of development branch dev-02 back to main trunk on 2005-10-25

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26