1 /*
2     Copyright (C) 1998  Dennis Roddeman
3     email: dennis.roddeman@feat.nl
4 
5     This program is free software; you can redistribute it and/or modify
6     it under the terms of the GNU General Public License as published by
7     the Free Software Foundation; either version 2 of the License, or
8     (at your option) any later version.
9 
10     This program is distributed in the hope that it will be useful,
11     but WITHOUT ANY WARRANTY; without even the implied warranty of
12     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
13     GNU General Public License for more details.
14 
15     You should have received a copy of the GNU General Public License
16     along with this program; if not, write to the Free Software
17     59 Temple Place, Suite 330, Boston, MA, 02111-1307, USA
18 */
19 
20 #include "tochnog.h"
21 
22 #define MNOD_MACRO 27
23 #define MEL_MACRO 8
24 #define EPSPHI 1.e-5
25 
macro(void)26 void macro( void )
27 
28 {
29   long int i=0, n=0, nthick=0, ncirc=0, ny=0, nz=0, nx=0, icontrol=0,
30     ix=0, ithick=0, icirc=0, iy=0, iz=0, inol=0, inod=0, macro_item=0,
31     element=0, idim=0,
32     length_mesh_macro=0, length_mesh_macro_parameters=0, isnot360=0,
33     use_control_refine_globally_geometry=0, nnod_macro=0, nel_macro=0,
34     length_control_refine_globally=0, max_geometry=0,ncircx=0,
35     node_boundary=-YES, length=0, nnol=0, element_group=0,
36     ok=0, geometry_macro_name=0, swit=0,
37     test1=0, test2=0, test3=0,
38     control_mesh_macro_element=0,
39     control_mesh_macro_set_node_boundary=-YES,
40     ldum=0, idum[1], node_remesh_allowed[MDIM], nodes[MNOL], el[1+MNOL],
41     control_refine_globally[4], control_refine_globally_geometry[2],
42     nodes_macro[MNOL][MEL_MACRO], *control_mesh_macro=NULL;
43   double radius=0., thickness=0., lenx=0., leny=0., lenz=0.,
44     r=0., phi=0., dx=0., dy=0., dz=0., phi0=0., phi1=2.0*PIRAD,
45     xuij=0., yuij=0., zuij=0., ij1=0., alpha=0., beta=0., ONE=0., TWO=0., THREE=0.,
46     ddum[1], work[MDIM], coord[MDIM], coordloc[MDIM], rot[MDIM*MDIM],
47     xc[MDIM], yc[MDIM], zc[MDIM], xc0[MDIM], yc0[MDIM], zc0[MDIM],
48     node_dof[MUKNWN], coord_macro[MDIM][MNOD_MACRO],
49     geometry_macro[MDIM+2], *control_mesh_macro_parameters=NULL;
50 
51   db( ICONTROL, 0, &icontrol, ddum, ldum, VERSION_NORMAL, GET );
52   if ( db_active_index( CONTROL_MESH_MACRO, icontrol, VERSION_NORMAL ) ) {
53     swit = set_swit(-1,-1,"macro");
54     if ( swit ) pri( "In routine MACRO" );
55     db( CONTROL_MESH_MACRO_SET_NODE_BOUNDARY, icontrol,
56       &control_mesh_macro_set_node_boundary, ddum, ldum, VERSION_NORMAL,
57       GET_IF_EXISTS );
58     control_mesh_macro =
59       db_int( CONTROL_MESH_MACRO, icontrol, VERSION_NORMAL );
60     length_mesh_macro = db_len( CONTROL_MESH_MACRO, icontrol, VERSION_NORMAL );
61     length_mesh_macro_parameters = db_len( CONTROL_MESH_MACRO_PARAMETERS,
62       icontrol, VERSION_NORMAL );
63     macro_item = control_mesh_macro[0];
64     element_group = control_mesh_macro[1];
65 
66     ok = 0;
67     if      ( macro_item==-BAR && length_mesh_macro==3 ) ok = 1;
68     else if ( macro_item==-RECTANGLE && length_mesh_macro==4 ) ok = 1;
69     else if ( macro_item==-CIRCLE && length_mesh_macro==3 ) ok = 1;
70     else if ( macro_item==-CIRCLE_HOLLOW && length_mesh_macro==4 ) ok = 1;
71     else if ( macro_item==-BRICK && length_mesh_macro==5 ) ok = 1;
72     else if ( macro_item==-CYLINDER_HOLLOW && length_mesh_macro==5 ) ok = 1;
73     else if ( macro_item==-SPHERE && length_mesh_macro==3 ) ok = 1;
74     if ( !ok ) db_error( CONTROL_MESH_MACRO, icontrol );
75 
76     ok = 0;
77     if      ( macro_item==-BAR && length_mesh_macro_parameters==2 ) ok = 1;
78     else if ( macro_item==-RECTANGLE && length_mesh_macro_parameters==4 ) ok = 1;
79     else if ( macro_item==-CIRCLE && length_mesh_macro_parameters==3 ) ok = 1;
80     else if ( macro_item==-CIRCLE_HOLLOW && length_mesh_macro_parameters==6 ) ok = 1;
81     else if ( macro_item==-BRICK && length_mesh_macro_parameters==6 ) ok = 1;
82     else if ( macro_item==-CYLINDER_HOLLOW && length_mesh_macro_parameters==10 ) ok = 1;
83     else if ( macro_item==-SPHERE && length_mesh_macro_parameters==4 ) ok = 1;
84     if ( !ok ) db_error( CONTROL_MESH_MACRO_PARAMETERS, icontrol );
85 
86     control_mesh_macro_parameters =
87       db_dbl( CONTROL_MESH_MACRO_PARAMETERS, icontrol, VERSION_NORMAL );
88     if ( macro_item==-CYLINDER_HOLLOW ) {
89       lenx = array_distance( &control_mesh_macro_parameters[0],
90         &control_mesh_macro_parameters[ndim], work, ndim );
91       xc[0] = control_mesh_macro_parameters[0];
92       yc[0] = control_mesh_macro_parameters[1];
93       zc[0] = control_mesh_macro_parameters[2];
94       xc[1] = control_mesh_macro_parameters[3];
95       yc[1] = control_mesh_macro_parameters[4];
96       zc[1] = control_mesh_macro_parameters[5];
97       array_set( xc0, 0., 2 );
98       array_set( yc0, 0., 2 );
99       array_set( zc0, 0., 2 );
100       xc0[1] = lenx;
101       xuij = (xc[1]-xc0[1]) - (xc[0]-xc0[0]);
102       yuij = (yc[1]-yc0[1]) - (yc[0]-yc0[0]);
103       zuij = (zc[1]-zc0[1]) - (zc[0]-zc0[0]);
104       ij1 = sqrt( scalar_square(lenx+xuij) + scalar_square(zuij) );
105       if ( scalar_dabs(ij1)<1.e-8 )
106         alpha = PIRAD/2.;
107       else
108         alpha = acos( (lenx+xuij)/ij1 );
109       beta = asin( yuij/lenx );
110       rot[0] = cos(alpha)*cos(beta);
111       rot[1] = sin(beta);
112       rot[2] = sin(alpha)*cos(beta);
113       rot[3] = -cos(alpha)*sin(beta);
114       rot[4] = cos(beta);
115       rot[5] = -sin(alpha)*sin(beta);
116       rot[6] = -sin(alpha);
117       rot[7] = 0.;
118       rot[8] = cos(alpha);
119     }
120     if      ( macro_item==-CIRCLE_HOLLOW || macro_item==-CYLINDER_HOLLOW ) {
121       if ( macro_item==-CIRCLE_HOLLOW ) {
122         nx = 1;
123         phi0=control_mesh_macro_parameters[4]*PIRAD/180.0;
124         phi1=control_mesh_macro_parameters[5]*PIRAD/180.0;
125         nthick = control_mesh_macro[2];
126         if ( nthick<2 ) db_error( CONTROL_MESH_MACRO, icontrol );
127         ncirc = control_mesh_macro[3];
128         if ( ncirc<2 ) db_error( CONTROL_MESH_MACRO, icontrol );
129         if ( scalar_dabs(phi1-phi0-2.0*PIRAD)<EPSPHI )
130           isnot360 = 0;
131         else
132           isnot360 = 1;
133         ncircx = ncirc + isnot360;
134         lenx = 0.;
135         radius = control_mesh_macro_parameters[2];
136         thickness = control_mesh_macro_parameters[3];
137       }
138       else {
139         assert( macro_item==-CYLINDER_HOLLOW );
140         nx = control_mesh_macro[2];
141         if ( nx<2 ) db_error( CONTROL_MESH_MACRO, icontrol );
142         nthick = control_mesh_macro[3];
143         if ( nthick<2 ) db_error( CONTROL_MESH_MACRO, icontrol );
144         ncirc = control_mesh_macro[4];
145         if ( ncirc<2 ) db_error( CONTROL_MESH_MACRO, icontrol );
146         radius = control_mesh_macro_parameters[6];
147         thickness = control_mesh_macro_parameters[7];
148         phi0=control_mesh_macro_parameters[8]*PIRAD/180.0;
149         phi1=control_mesh_macro_parameters[9]*PIRAD/180.0;
150         if ( scalar_dabs(phi1-phi0-2.0*PIRAD)<EPSPHI )
151           isnot360 = 0;
152         else
153           isnot360 = 1;
154         ncircx = ncirc + isnot360;
155       }
156       for ( ix=0; ix<nx; ix++ ) {
157         dx = ix*lenx/(nx-1);
158         for ( icirc=0; icirc<ncircx; icirc++ ) {
159           phi = icirc*(phi1-phi0)/(ncirc)+phi0;
160           for ( ithick=0; ithick<nthick; ithick++ ) {
161             r = radius - thickness/2. + ithick*thickness/(nthick-1);
162             if      ( ndim==2 ) {
163               coordloc[0] = r*cos(phi);
164               coordloc[1] = r*sin(phi);
165               array_move( coordloc, coord, ndim );
166             }
167             else {
168               assert( ndim==3 );
169               coordloc[0] = dx;
170               coordloc[1] = r*cos(phi);
171               coordloc[2] = r*sin(phi);
172               matrix_atb( rot, coordloc, coord, ndim, ndim, 1 );
173             }
174             array_add( coord, control_mesh_macro_parameters, coord, ndim );
175             db( NODE, inod, idum, coord, ndim, VERSION_MACRO, PUT );
176             db( NODE_START_REFINED, inod, idum, coord, ndim, VERSION_MACRO, PUT );
177             length = 1; db( NODE_MACRO_GENERATE, inod, &icontrol, ddum, length, VERSION_MACRO, PUT );
178             if ( nuknwn>0 ) {
179               array_set( node_dof, 0., nuknwn );
180               db( NODE_DOF, inod, idum, node_dof, nuknwn, VERSION_MACRO, PUT );
181               db( NODE_DOF_START_REFINED, inod, idum, node_dof, nuknwn, VERSION_MACRO, PUT );
182             }
183             if ( ithick==0 || ithick==nthick-1 || (isnot360==1 && (icirc==0 || icirc==ncirc)) ||
184                 ( ndim==3 && (ix==0 || ix==nx-1) ) ) {
185               if (  control_mesh_macro_set_node_boundary==-YES) {
186                 length = 1; db( NODE_BOUNDARY, inod, &node_boundary, ddum,
187                   length, VERSION_MACRO, PUT );
188               }
189             }
190             inod++;
191           }
192         }
193       }
194       for ( ix=0; ix<scalar_imax(1,nx-1); ix++ ) {
195         for ( icirc=0; icirc<ncirc; icirc++ ) {
196           for ( ithick=0; ithick<nthick-1; ithick++ ) {
197             if ( ndim==2 ) {
198               el[0] = -QUAD4;
199               nnol = 4;
200               nodes[0] = icirc*nthick+ithick;
201               nodes[1] = icirc*nthick+ithick + 1;
202               if ( icirc==ncirc-1 && isnot360==0 ) {
203                 nodes[2] = 0*nthick+ithick;
204                 nodes[3] = 0*nthick+ithick + 1;
205               }
206               else {
207                 nodes[2] = (icirc+1)*nthick+ithick;
208                 nodes[3] = (icirc+1)*nthick+ithick + 1;
209               }
210             }
211             else {
212               assert( ndim==3 );
213               el[0] = -HEX8;
214               nnol = 8;
215               nodes[0] = ix*ncircx*nthick+icirc*nthick+ithick;
216               nodes[1] = ix*ncircx*nthick+icirc*nthick+ithick + 1;
217               nodes[4] = (ix+1)*ncircx*nthick+icirc*nthick+ithick;
218               nodes[5] = (ix+1)*ncircx*nthick+icirc*nthick+ithick + 1;
219               if ( icirc==ncirc-1 && isnot360==0 ) {
220                 nodes[2] = ix*ncircx*nthick+0*nthick+ithick;
221                 nodes[3] = ix*ncircx*nthick+0*nthick+ithick + 1;
222                 nodes[6] = (ix+1)*ncircx*nthick+0*nthick+ithick;
223                 nodes[7] = (ix+1)*ncircx*nthick+0*nthick+ithick + 1;
224               }
225               else {
226                 nodes[2] = ix*ncircx*nthick+(icirc+1)*nthick+ithick;
227                 nodes[3] = ix*ncircx*nthick+(icirc+1)*nthick+ithick + 1;
228                 nodes[6] = (ix+1)*ncircx*nthick+(icirc+1)*nthick+ithick;
229                 nodes[7] = (ix+1)*ncircx*nthick+(icirc+1)*nthick+ithick + 1;
230               }
231             }
232 
233             array_move( nodes, &el[1], nnol );
234             length = 1 + nnol;
235             db( ELEMENT, element, el, ddum, length, VERSION_MACRO, PUT );
236             length = 1;
237             db( ELEMENT_GROUP, element, &element_group, ddum, length, VERSION_MACRO, PUT );
238             length = 1;
239             db( ELEMENT_MACRO_GENERATE, element, &icontrol, ddum, length, VERSION_MACRO, PUT );
240             element++;
241           }
242         }
243       }
244     }
245     else if ( macro_item==-SPHERE || macro_item==-CIRCLE ) {
246       ONE = 1.;
247       TWO = 1./sqrt(2.);
248       THREE = 1./sqrt(3.);
249       n = control_mesh_macro[2];
250       if ( n<0 ) db_error( CONTROL_MESH_MACRO, icontrol );
251       radius = control_mesh_macro_parameters[ndim];
252       array_set( &coord_macro[0][0], 0., MDIM*MNOD_MACRO );
253       if ( macro_item==-SPHERE ) {
254         nnod_macro = 27;
255         nel_macro = 8;
256         coord_macro[0][0]   = -THREE; coord_macro[1][0]   = -THREE; coord_macro[2][0]   = -THREE;
257         coord_macro[0][1]   =  0.;    coord_macro[1][1]   = -TWO;   coord_macro[2][1]   = -TWO;
258         coord_macro[0][2]   = +THREE; coord_macro[1][2]   = -THREE; coord_macro[2][2]   = -THREE;
259         coord_macro[0][3]   = -TWO;   coord_macro[1][3]   =  0.;    coord_macro[2][3]   = -TWO;
260         coord_macro[0][4]   =  0.;    coord_macro[1][4]   = -0.;    coord_macro[2][4]   = -ONE;
261         coord_macro[0][5]   = +TWO;   coord_macro[1][5]   =  0.;    coord_macro[2][5]   = -TWO;
262         coord_macro[0][6]   = -THREE; coord_macro[1][6]   = +THREE; coord_macro[2][6]   = -THREE;
263         coord_macro[0][7]   =  0.;    coord_macro[1][7]   = +TWO;   coord_macro[2][7]   = -TWO;
264         coord_macro[0][8]   = +THREE; coord_macro[1][8]   = +THREE; coord_macro[2][8]   = -THREE;
265 
266         coord_macro[0][9]   = -TWO;   coord_macro[1][9]   = -TWO;   coord_macro[2][9]   =  0.;
267         coord_macro[0][10]  =  0.;    coord_macro[1][10]  = -ONE;   coord_macro[2][10]  =  0.;
268         coord_macro[0][11]  = +TWO;   coord_macro[1][11]  = -TWO;   coord_macro[2][11]  =  0.;
269         coord_macro[0][12]  = -ONE;   coord_macro[1][12]  =  0.;    coord_macro[2][12]  =  0.;
270         coord_macro[0][13]  =  0.;    coord_macro[1][13]  =  0.;    coord_macro[2][13]  =  0.;
271         coord_macro[0][14]  = +ONE;   coord_macro[1][14]  =  0.;    coord_macro[2][14]  =  0.;
272         coord_macro[0][15]  = -TWO;   coord_macro[1][15]  = +TWO;   coord_macro[2][15]  =  0.;
273         coord_macro[0][16]  =  0.;    coord_macro[1][16]  = +ONE;   coord_macro[2][16]  =  0.;
274         coord_macro[0][17]  = +TWO;   coord_macro[1][17]  = +TWO;   coord_macro[2][17]  =  0.;
275 
276         coord_macro[0][18]  = -THREE; coord_macro[1][18]  = -THREE; coord_macro[2][18]  = +THREE;
277         coord_macro[0][19]  =  0.;    coord_macro[1][19]  = -TWO;   coord_macro[2][19]  = +TWO;
278         coord_macro[0][20]  = +THREE; coord_macro[1][20]  = -THREE; coord_macro[2][20]  = +THREE;
279         coord_macro[0][21]  = -TWO;   coord_macro[1][21]  =  0.;    coord_macro[2][21]  = +TWO;
280         coord_macro[0][22]  =  0.;    coord_macro[1][22]  =  0.;    coord_macro[2][22]  = +ONE;
281         coord_macro[0][23]  = +TWO;   coord_macro[1][23]  =  0.;    coord_macro[2][23]  = +TWO;
282         coord_macro[0][24]  = -THREE; coord_macro[1][24]  = +THREE; coord_macro[2][24]  = +THREE;
283         coord_macro[0][25]  =  0.;    coord_macro[1][25]  = +TWO;   coord_macro[2][25]  = +TWO;
284         coord_macro[0][26]  = +THREE; coord_macro[1][26]  = +THREE; coord_macro[2][26]  = +THREE;
285 
286         nodes_macro[0][0]  = 0;  nodes_macro[1][0]  = 3;  nodes_macro[2][0]  = 1;  nodes_macro[3][0]  = 4;
287         nodes_macro[4][0]  = 9;  nodes_macro[5][0]  = 12; nodes_macro[6][0]  = 10; nodes_macro[7][0]  = 13;
288 
289         nodes_macro[0][1]  = 1;  nodes_macro[1][1]  = 4;  nodes_macro[2][1]  = 2;  nodes_macro[3][1]  = 5;
290         nodes_macro[4][1]  = 10; nodes_macro[5][1]  = 13; nodes_macro[6][1]  = 11; nodes_macro[7][1]  = 14;
291 
292         nodes_macro[0][2]  = 3;  nodes_macro[1][2]  = 6;  nodes_macro[2][2]  = 4;  nodes_macro[3][2]  = 7;
293         nodes_macro[4][2]  = 12; nodes_macro[5][2]  = 15; nodes_macro[6][2]  = 13; nodes_macro[7][2]  = 16;
294 
295         nodes_macro[0][3]  = 4;  nodes_macro[1][3]  = 7;  nodes_macro[2][3]  = 5;  nodes_macro[3][3]  = 8;
296         nodes_macro[4][3]  = 13; nodes_macro[5][3]  = 16; nodes_macro[6][3]  = 14; nodes_macro[7][3]  = 17;
297 
298         nodes_macro[0][4]  = 9;  nodes_macro[1][4]  = 12; nodes_macro[2][4]  = 10; nodes_macro[3][4]  = 13;
299         nodes_macro[4][4]  = 18; nodes_macro[5][4]  = 21; nodes_macro[6][4]  = 19; nodes_macro[7][4]  = 22;
300 
301         nodes_macro[0][5]  = 10; nodes_macro[1][5]  = 13; nodes_macro[2][5]  = 11; nodes_macro[3][5]  = 14;
302         nodes_macro[4][5]  = 19; nodes_macro[5][5]  = 22; nodes_macro[6][5]  = 20; nodes_macro[7][5]  = 23;
303 
304         nodes_macro[0][6]  = 12; nodes_macro[1][6]  = 15; nodes_macro[2][6]  = 13; nodes_macro[3][6]  = 16;
305         nodes_macro[4][6]  = 21; nodes_macro[5][6]  = 24; nodes_macro[6][6]  = 22; nodes_macro[7][6]  = 25;
306 
307         nodes_macro[0][7]  = 13; nodes_macro[1][7]  = 16; nodes_macro[2][7]  = 14; nodes_macro[3][7]  = 17;
308         nodes_macro[4][7]  = 22; nodes_macro[5][7]  = 25; nodes_macro[6][7]  = 23; nodes_macro[7][7]  = 26;
309 
310       }
311       else {
312         assert( macro_item==-CIRCLE );
313         nnod_macro = 9;
314         nel_macro = 4;
315 
316         coord_macro[0][0]   = -TWO; coord_macro[1][0]   = -TWO;
317         coord_macro[0][1]   =   0.; coord_macro[1][1]   = -ONE;
318         coord_macro[0][2]   = +TWO; coord_macro[1][2]   = -TWO;
319 
320         coord_macro[0][3]   = -ONE; coord_macro[1][3]   =    0.;
321         coord_macro[0][4]   =   0.; coord_macro[1][4]   =    0.;
322         coord_macro[0][5]   = +ONE; coord_macro[1][5]   =    0.;
323 
324         coord_macro[0][6]   = -TWO; coord_macro[1][6]   = +TWO;
325         coord_macro[0][7]   =   0.; coord_macro[1][7]   = +ONE;
326         coord_macro[0][8]   = +TWO; coord_macro[1][8]   = +TWO;
327 
328         nodes_macro[0][0]  = 0;  nodes_macro[1][0]  = 1;  nodes_macro[2][0]  = 3;  nodes_macro[3][0]  = 4;
329 
330         nodes_macro[0][1]  = 1;  nodes_macro[1][1]  = 2;  nodes_macro[2][1]  = 4;  nodes_macro[3][1]  = 5;
331 
332         nodes_macro[0][2]  = 3;  nodes_macro[1][2]  = 4;  nodes_macro[2][2]  = 6;  nodes_macro[3][2]  = 7;
333 
334         nodes_macro[0][3]  = 4;  nodes_macro[1][3]  = 5;  nodes_macro[2][3]  = 7;  nodes_macro[3][3]  = 8;
335       }
336       array_multiply( &coord_macro[0][0], &coord_macro[0][0], radius, MDIM*MNOD_MACRO );
337 
338       for ( inod=0; inod<nnod_macro; inod++ ) {
339         for ( idim=0; idim<ndim; idim++ ) coordloc[idim] = coord_macro[idim][inod];
340         array_move( coordloc, coord, ndim );
341         array_add( coord, control_mesh_macro_parameters, coord, ndim );
342         db( NODE, inod, idum, coord, ndim, VERSION_MACRO, PUT );
343         db( NODE_START_REFINED, inod, idum, coord, ndim, VERSION_MACRO, PUT );
344         length = 1; db( NODE_MACRO_GENERATE, inod, &icontrol, ddum, length, VERSION_MACRO, PUT );
345         if ( nuknwn>0 ) {
346           array_set( node_dof, 0., nuknwn );
347           db( NODE_DOF, inod, idum, node_dof, nuknwn, VERSION_MACRO, PUT );
348           db( NODE_DOF_START_REFINED, inod, idum, node_dof, nuknwn, VERSION_MACRO, PUT );
349         }
350         if ( control_mesh_macro_set_node_boundary==-YES && inod!=(nnod_macro-1)/2 ) {
351           length = 1; db( NODE_BOUNDARY, inod, &node_boundary, ddum,
352             length, VERSION_MACRO, PUT );
353         }
354       }
355       for ( element=0; element<nel_macro; element++ ) {
356         if ( ndim==3 ) {
357           el[0] = -HEX8; nnol = 8;
358           geometry_macro_name = -GEOMETRY_SPHERE;
359         }
360         else {
361           assert( ndim==2 );
362           el[0] = -QUAD4; nnol = 4;
363           geometry_macro_name = -GEOMETRY_CIRCLE;
364         }
365         for ( inol=0; inol<nnol; inol++ ) el[1+inol] = nodes_macro[inol][element];
366         length = 1 + nnol;
367         db( ELEMENT, element, el, ddum, length, VERSION_MACRO, PUT );
368         length = 1;
369         db( ELEMENT_GROUP, element, &element_group, ddum, length, VERSION_MACRO, PUT );
370         length = 1;
371         db( ELEMENT_MACRO_GENERATE, element, &icontrol, ddum, length, VERSION_MACRO, PUT );
372       }
373 
374       array_move( control_mesh_macro_parameters, geometry_macro, ndim+1 );
375       geometry_macro[ndim+1] = EPS_COORD;
376       db_max_index( geometry_macro_name, max_geometry, VERSION_NORMAL, GET );
377       max_geometry++; length = ndim + 2;
378       db( geometry_macro_name, max_geometry, idum, geometry_macro, length, VERSION_NORMAL, PUT );
379       use_control_refine_globally_geometry = 1;
380       control_refine_globally_geometry[0] = geometry_macro_name;
381       control_refine_globally_geometry[1] = max_geometry;
382       for ( i=0; i<n; i++ ) {
383         length_control_refine_globally = 1;
384         control_refine_globally[0] = -H_REFINEMENT;
385         refine_globally( control_refine_globally, length_control_refine_globally,
386           use_control_refine_globally_geometry, control_refine_globally_geometry,
387           PROJECT_EXACT, VERSION_MACRO );
388       }
389       db_delete_index( geometry_macro_name, max_geometry, VERSION_NORMAL );
390 
391     }
392     else if ( macro_item==-BAR || macro_item==-RECTANGLE || macro_item==-BRICK ) {
393       if      ( macro_item==-BAR ) {
394         nx = 1;
395         ny = 1;
396         nz = control_mesh_macro[2];
397         if ( nz<2 ) db_error( CONTROL_MESH_MACRO, icontrol );
398         lenz = control_mesh_macro_parameters[ndim];
399       }
400       else if ( macro_item==-RECTANGLE ) {
401         nx = 1;
402         ny = control_mesh_macro[2];
403         if ( ny<2 ) db_error( CONTROL_MESH_MACRO, icontrol );
404         nz = control_mesh_macro[3];
405         if ( nz<2 ) db_error( CONTROL_MESH_MACRO, icontrol );
406         leny = control_mesh_macro_parameters[ndim];
407         lenz = control_mesh_macro_parameters[ndim+1];
408       }
409       else {
410         assert( macro_item==-BRICK );
411         nx = control_mesh_macro[2];
412         if ( nx<2 ) db_error( CONTROL_MESH_MACRO, icontrol );
413         ny = control_mesh_macro[3];
414         if ( ny<2 ) db_error( CONTROL_MESH_MACRO, icontrol );
415         nz = control_mesh_macro[4];
416         if ( nz<2 ) db_error( CONTROL_MESH_MACRO, icontrol );
417         lenx = control_mesh_macro_parameters[ndim];
418         leny = control_mesh_macro_parameters[ndim+1];
419         lenz = control_mesh_macro_parameters[ndim+2];
420       }
421       for ( ix=0; ix<nx; ix++ ) {
422         if ( ndim>=3 ) dx = ix*lenx/(nx-1);
423         for ( iy=0; iy<ny; iy++ ) {
424           if ( ndim>=2 ) dy = iy*leny/(ny-1);
425           for ( iz=0; iz<nz; iz++ ) {
426             if ( ndim>=1 ) dz = iz*lenz/(nz-1);
427             if      ( ndim==1 ) {
428               coord[0] = -0.5*lenz + dz;
429             }
430             else if ( ndim==2 ) {
431               coord[0] = -0.5*leny + dy;
432               coord[1] = -0.5*lenz + dz;
433             }
434             else {
435               assert( ndim==3 );
436               coord[0] = -0.5*lenx + dx;
437               coord[1] = -0.5*leny + dy;
438               coord[2] = -0.5*lenz + dz;
439             }
440             array_add( coord, control_mesh_macro_parameters, coord, ndim );
441             db( NODE, inod, idum, coord, ndim, VERSION_MACRO, PUT );
442             db( NODE_START_REFINED, inod, idum, coord, ndim, VERSION_MACRO, PUT );
443             length = 1; db( NODE_MACRO_GENERATE, inod, &icontrol, ddum, length, VERSION_MACRO, PUT );
444             if ( nuknwn>0 ) {
445               array_set( node_dof, 0., nuknwn );
446               db( NODE_DOF, inod, idum, node_dof, nuknwn, VERSION_MACRO, PUT );
447               db( NODE_DOF_START_REFINED, inod, idum, node_dof, nuknwn, VERSION_MACRO, PUT );
448             }
449             test1 = ndim>=3 && (ix==0 || ix==nx-1);
450             test2 = ndim>=2 && (iy==0 || iy==ny-1);
451             test3 = ndim>=1 && (iz==0 || iz==nz-1);
452             if ( test1 || test2 || test3 ) {
453               if (  control_mesh_macro_set_node_boundary==-YES ) {
454                 length = 1; db( NODE_BOUNDARY, inod, &node_boundary, ddum,
455                   length, VERSION_MACRO, PUT );
456               }
457             }
458             array_set( node_remesh_allowed, -YES, ndim );
459             if ( test1 ) node_remesh_allowed[0] = -NO;
460             if ( test2 ) node_remesh_allowed[1] = -NO;
461             if ( test3 ) node_remesh_allowed[2] = -NO;
462             length = ndim; db( NODE_REMESH_ALLOWED, inod, node_remesh_allowed, ddum,
463                length, VERSION_MACRO, PUT );
464             inod++;
465           }
466         }
467       }
468       for ( ix=0; ix<scalar_imax(1,nx-1); ix++ ) {
469         for ( iy=0; iy<scalar_imax(1,ny-1); iy++ ) {
470           for ( iz=0; iz<scalar_imax(1,nz-1); iz++ ) {
471             if      ( ndim==1 ) {
472               el[0] = -BAR2;
473               nnol = 2;
474               nodes[0] = iz;
475               nodes[1] = iz + 1;
476             }
477             else if ( ndim==2 ) {
478               el[0] = -QUAD4;
479               nnol = 4;
480               nodes[0] = iy*nz+iz;
481               nodes[1] = iy*nz+iz + 1;
482               nodes[2] = (iy+1)*nz+iz;
483               nodes[3] = (iy+1)*nz+iz + 1;
484             }
485             else {
486               assert( ndim==3 );
487               el[0] = -HEX8;
488               nnol = 8;
489               nodes[0] = ix*nz*ny+iy*nz+iz;
490               nodes[1] = ix*nz*ny+iy*nz+iz + 1;
491               nodes[2] = ix*nz*ny+(iy+1)*nz+iz;
492               nodes[3] = ix*nz*ny+(iy+1)*nz+iz + 1;
493               nodes[4] = (ix+1)*nz*ny+iy*nz+iz;
494               nodes[5] = (ix+1)*nz*ny+iy*nz+iz + 1;
495               nodes[6] = (ix+1)*nz*ny+(iy+1)*nz+iz;
496               nodes[7] = (ix+1)*nz*ny+(iy+1)*nz+iz + 1;
497             }
498             array_move( nodes, &el[1], nnol );
499             length = 1 + nnol;
500             db( ELEMENT, element, el, ddum, length, VERSION_MACRO, PUT );
501             length = 1;
502             db( ELEMENT_GROUP, element, &element_group, ddum, length, VERSION_MACRO, PUT );
503             length = 1;
504             db( ELEMENT_MACRO_GENERATE, element, &icontrol, ddum, length, VERSION_MACRO, PUT );
505             element++;
506           }
507         }
508       }
509     }
510     else
511       db_error( CONTROL_MESH_MACRO, icontrol );
512 
513       // split or p_refine mesh
514     if ( db( CONTROL_MESH_MACRO_ELEMENT, icontrol, &control_mesh_macro_element,
515         ddum, ldum, VERSION_NORMAL, GET_IF_EXISTS ) ) {
516       if      ( control_mesh_macro_element==-TRIA3 )
517         mesh_split( VERSION_MACRO );
518       else if ( control_mesh_macro_element==-TET4 )
519         mesh_split( VERSION_MACRO );
520       else if ( control_mesh_macro_element==-BAR3  ||
521                 control_mesh_macro_element==-QUAD9 ||
522                 control_mesh_macro_element==-HEX27 ) {
523         use_control_refine_globally_geometry = 0;
524         length_control_refine_globally = 1;
525         control_refine_globally[0] = -P_REFINEMENT;
526         refine_globally( control_refine_globally, length_control_refine_globally,
527           use_control_refine_globally_geometry, control_refine_globally_geometry,
528           PROJECT_EXACT, VERSION_MACRO );
529       }
530       else if ( control_mesh_macro_element==-TRIA6  ||
531                 control_mesh_macro_element==-TET10 ) {
532         use_control_refine_globally_geometry = 0;
533         length_control_refine_globally = 1;
534         control_refine_globally[0] = -P_REFINEMENT;
535         refine_globally( control_refine_globally, length_control_refine_globally,
536           use_control_refine_globally_geometry, control_refine_globally_geometry,
537           PROJECT_EXACT, VERSION_MACRO );
538         mesh_split( VERSION_MACRO );
539       }
540     }
541 
542       // renumber
543     renumbering( VERSION_MACRO, NO, 1, 1, idum, idum );
544 
545       // add to existing mesh
546     mesh_add( VERSION_MACRO, VERSION_NORMAL );
547 
548     if ( swit ) pri( "Out routine MACRO" );
549 
550   }
551 
552 }
553 
554