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

Annotation of /trunk/esys2/finley/src/finleyC/Mesh_rec4.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 149 - (hide annotations)
Thu Sep 1 03:31:39 2005 UTC (14 years, 3 months ago) by jgs
File MIME type: text/plain
File size: 8441 byte(s)
Merge of development branch dev-02 back to main trunk on 2005-09-01

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26