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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 153 - (show annotations)
Tue Oct 25 01:51:20 2005 UTC (13 years, 11 months 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 /*
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,M0,M1;
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 if (N0<=N1) {
46 M0=1;
47 M1=N0;
48 } else {
49 M0=N1;
50 M1=1;
51 }
52
53 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 if (! Finley_noError()) return NULL;
72
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 if (! Finley_noError()) {
83 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 if (! Finley_noError()) {
93 Finley_Mesh_dealloc(out);
94 return NULL;
95 }
96 /* set nodes: */
97
98 #pragma omp parallel for private(i0,i1,k)
99 for (i1=0;i1<N1;i1++) {
100 for (i0=0;i0<N0;i0++) {
101 k=M0*i0+M1*i1;
102 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 out->Nodes->Id[k]=i0+N0*i1;
105 out->Nodes->Tag[k]=0;
106 out->Nodes->degreeOfFreedom[k]=M0*(i0%NDOF0) +M1*(i1%NDOF1);
107 }
108 }
109 /* tags for the faces: */
110 if (!periodic[1]) {
111 for (i0=0;i0<N0;i0++) {
112 out->Nodes->Tag[M0*i0+M1*0]+=10;
113 out->Nodes->Tag[M0*i0+M1*(N1-1)]+=20;
114 }
115 }
116 if (!periodic[0]) {
117 for (i1=0;i1<N1;i1++) {
118 out->Nodes->Tag[M0*0+M1*i1]+=1;
119 out->Nodes->Tag[M0*(N0-1)+M1*i1]+=2;
120 }
121 }
122
123 /* 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 out->Elements->minColor=0;
143 out->Elements->maxColor=COLOR_MOD(0)+3*COLOR_MOD(0);
144
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 out->FaceElements->minColor=0;
253 out->FaceElements->maxColor=7;
254
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 #ifdef Finley_TRACE
266 printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);
267 #endif
268
269 if (! Finley_noError()) {
270 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