1 /*
2 * ***************************************************************************
3 * MC = < Manifold Code >
4 * Copyright (C) 1994-- Michael Holst
5 *
6 * This library is free software; you can redistribute it and/or
7 * modify it under the terms of the GNU Lesser General Public
8 * License as published by the Free Software Foundation; either
9 * version 2.1 of the License, or (at your option) any later version.
10 *
11 * This library is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14 * Lesser General Public License for more details.
15 *
16 * You should have received a copy of the GNU Lesser General Public
17 * License along with this library; if not, write to the Free Software
18 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19 *
20 * rcsid="$Id: aprx.c,v 1.29 2010/08/12 05:18:15 fetk Exp $"
21 * ***************************************************************************
22 */
23
24 /*
25 * ***************************************************************************
26 * File: aprx.c
27 *
28 * Purpose: Class Aprx: methods.
29 *
30 * Author: Michael Holst
31 * ***************************************************************************
32 */
33
34 #include "aprx_p.h"
35
36 VEMBED(rcsid="$Id: aprx.c,v 1.29 2010/08/12 05:18:15 fetk Exp $")
37
38 /*
39 * ***************************************************************************
40 * Class Aprx: Inlineable methods
41 * ***************************************************************************
42 */
43 #if !defined(VINLINE_APRX)
44 #endif /* if !defined(VINLINE_GEM) */
45
46 /*
47 * ***************************************************************************
48 * Class Aprx: Non-inlineable methods
49 * ***************************************************************************
50 */
51
52 /*
53 * ***************************************************************************
54 * Routine: Aprx_ctor
55 *
56 * Purpose: Construct a linear Approximation object
57 *
58 * Author: Michael Holst
59 * ***************************************************************************
60 */
Aprx_ctor(Vmem * vmem,Gem * tgm,PDE * tpde)61 VPUBLIC Aprx* Aprx_ctor(Vmem *vmem, Gem *tgm, PDE *tpde)
62 {
63 /* Create a linear approximation */
64 return Aprx_ctor2(vmem,tgm,tpde,1);
65 }
66
67 /*
68 * ***************************************************************************
69 * Routine: Aprx_ctor2
70 *
71 * Purpose: Approximation object constructor.
72 *
73 * Author: Michael Holst
74 * ***************************************************************************
75 */
Aprx_ctor2(Vmem * vmem,Gem * tgm,PDE * tpde,int order)76 VPUBLIC Aprx* Aprx_ctor2(Vmem *vmem, Gem *tgm, PDE *tpde,int order)
77 {
78 int i;
79 Aprx* thee = VNULL;
80
81 thee = Vmem_malloc( VNULL, 1, sizeof(Aprx) );
82 if (vmem == VNULL) {
83 thee->vmem = Vmem_ctor( "Aprx" );
84 thee->iMadeVmem = 1;
85 } else {
86 thee->vmem = vmem;
87 thee->iMadeVmem = 0;
88 }
89
90 VDEBUGIO("Aprx_ctor: CREATING object..");
91 VDEBUGIO("..done.\n");
92
93 thee->gm = tgm;
94 thee->pde = tpde;
95
96 thee->numB = PDE_vec(thee->pde);
97 thee->lnkG = Vset_ctor( thee->vmem, "G", sizeof( Link ), VMAX_OBJECTS, 0 );
98
99 thee->B = VNULL;
100 for (i=0; i<MAXV; i++) {
101 thee->re[i] = VNULL;
102 thee->reB[i] = VNULL;
103 }
104
105 thee->wev = VNULL;
106 thee->P = VNULL;
107
108 thee->order = order;
109
110 Aprx_reset(thee);
111
112 return thee;
113 }
114
115 /*
116 * ***************************************************************************
117 * Routine: Aprx_dtor
118 *
119 * Purpose: Approximation object destructor.
120 *
121 * Author: Michael Holst
122 * ***************************************************************************
123 */
Aprx_dtor(Aprx ** thee)124 VPUBLIC void Aprx_dtor(Aprx **thee)
125 {
126 int i;
127
128 VASSERT( (*thee) != VNULL );
129 if ((*thee) != VNULL) {
130
131 if ((*thee)->lnkG != VNULL) {
132 Vset_dtor( &((*thee)->lnkG) );
133 }
134
135 if ((*thee)->P != VNULL) {
136 Bmat_dtor( &((*thee)->P) );
137 }
138
139 if ((*thee)->B != VNULL) {
140 Bnode_dtor( &((*thee)->B) );
141 }
142 for (i=0; i<(*thee)->numB; i++) {
143 if ((*thee)->numBV[i] > 0) {
144 Vmem_free( ((*thee)->vmem), (*thee)->numBV[i], sizeof(int),
145 (void**)&((*thee)->ibv[i]) );
146 Vmem_free( ((*thee)->vmem), (*thee)->numBV[i], sizeof(double),
147 (void**)&((*thee)->bv[i]) );
148 }
149 }
150
151 for (i=0; i<MAXV; i++) {
152 if ((*thee)->re[i] != VNULL) {
153 Re_dtor( &((*thee)->re[i]) );
154 }
155 if ((*thee)->reB[i] != VNULL) {
156 Re_dtor( &((*thee)->reB[i]) );
157 }
158 }
159
160 if ((*thee)->wev != VNULL) {
161 Bvec_dtor( &((*thee)->wev) );
162 }
163
164 VDEBUGIO("Aprx_dtor: DESTROYING object..");
165 if ((*thee)->iMadeVmem) Vmem_dtor( &((*thee)->vmem) );
166 Vmem_free( VNULL, 1, sizeof(Aprx), (void**)thee );
167 VDEBUGIO("..done.\n");
168
169 (*thee) = VNULL;
170 }
171 }
172
173 /*
174 * ***************************************************************************
175 * Routine: Aprx_numB
176 *
177 * Purpose: Return the number of blocks.
178 *
179 * Author: Michael Holst
180 * ***************************************************************************
181 */
Aprx_numB(Aprx * thee)182 VPUBLIC int Aprx_numB(Aprx *thee)
183 {
184 return thee->numB;
185 }
186
187 /*
188 * ***************************************************************************
189 * Routine: Aprx_P
190 *
191 * Purpose: Return a pointer to the current prolongation matrix.
192 *
193 * Author: Michael Holst
194 * ***************************************************************************
195 */
Aprx_P(Aprx * thee)196 VPUBLIC Bmat *Aprx_P(Aprx *thee)
197 {
198 return thee->P;
199 }
200
201 /*
202 * ***************************************************************************
203 * Routine: Aprx_read
204 *
205 * Purpose: Read in the user-specified initial mesh given in the
206 * "MCSF" or "MCEF" format, and transform into our internal
207 * datastructures.
208 *
209 * Do a little more than a "Gem_read", in that we also
210 * initialize the extrinsic and intrinsic spatial dimensions
211 * corresponding to the input mesh, and we also then build the
212 * reference elements.
213 *
214 * Notes: See the documentation to Gem_read for a description of the
215 * mesh input data file format.
216 *
217 * Author: Michael Holst
218 * ***************************************************************************
219 */
Aprx_read(Aprx * thee,int key,Vio * sock)220 VPUBLIC int Aprx_read(Aprx *thee, int key, Vio *sock)
221 {
222 /* read in the mesh; return error if problems */
223 if ( !Gem_read(thee->gm,key,sock) ) {
224 return 0;
225 }
226
227 /* create node information for the mesh */
228 Aprx_reset(thee);
229 Aprx_buildNodes(thee);
230
231 /* error-free return */
232 return 1;
233 }
234
235 /*
236 * ***************************************************************************
237 * Routine: Aprx_reset
238 *
239 * Purpose: Reset all of the datastructures.
240 *
241 * Author: Michael Holst (updated by Ryan Szypowski)
242 * ***************************************************************************
243 */
Aprx_reset(Aprx * thee)244 VPUBLIC void Aprx_reset(Aprx *thee)
245 {
246 int i;
247
248 /* now reset some other things */
249 if ( thee->P != VNULL ) Bmat_dtor( &(thee->P) );
250 if ( thee->B != VNULL ) Bnode_dtor( &(thee->B) );
251 thee->P = VNULL;
252 thee->B = VNULL;
253
254 Vset_reset( thee->lnkG );
255
256 if (thee->I_V2E != VNULL) {
257 Vmem_free(thee->vmem,2*(thee->numE),sizeof(int),(void**)&(thee->I_V2E));
258 }
259 if (thee->I_V2F != VNULL) {
260 Vmem_free(thee->vmem,3*(thee->numF),sizeof(int),(void**)&(thee->I_V2F));
261 }
262 if (thee->I_V2S != VNULL) {
263 Vmem_free(thee->vmem,4*(thee->numS),sizeof(int),(void**)&(thee->I_V2S));
264 }
265 if (thee->I_E2F != VNULL) {
266 Vmem_free(thee->vmem,3*(thee->numF),sizeof(int),(void**)&(thee->I_E2F));
267 }
268 if (thee->I_E2S != VNULL) {
269 Vmem_free(thee->vmem,6*(thee->numS),sizeof(int),(void**)&(thee->I_E2S));
270 }
271 if (thee->I_F2S != VNULL) {
272 Vmem_free(thee->vmem,4*(thee->numS),sizeof(int),(void**)&(thee->I_F2S));
273 }
274 thee->numV = 0;
275 thee->numE = 0;
276 thee->numF = 0;
277 thee->numS = 0;
278
279 /* build new master elements for this mesh */
280 for (i=0; i<thee->numB; i++) {
281 if (thee->re[i] != VNULL) {
282 Re_dtor( &(thee->re[i]) );
283 }
284 thee->re[i] = Re_ctor(0, PDE_dim(thee->pde),
285 thee->pde->simplexBasisInit, thee->pde->simplexBasisForm,thee->order);
286 if (thee->reB[i] != VNULL) {
287 Re_dtor( &(thee->reB[i]) );
288 }
289 if (PDE_dim(thee->pde) == 3) {
290 /* Use cubic face bumps for 3d. */
291 thee->reB[i] = Re_ctor(1, PDE_dim(thee->pde),
292 thee->pde->simplexBasisInit, thee->pde->simplexBasisForm,3);
293 }
294 else {
295 thee->reB[i] = Re_ctor(0, PDE_dim(thee->pde),
296 thee->pde->simplexBasisInit, thee->pde->simplexBasisForm,thee->order);
297 }
298 }
299 }
300
301 /*
302 * ***************************************************************************
303 * Routine: Aprx_initNodes
304 *
305 * Purpose: Make node information.
306 *
307 * Author: Michael Holst
308 * ***************************************************************************
309 */
Aprx_initNodes(Aprx * thee)310 VPUBLIC void Aprx_initNodes(Aprx *thee)
311 {
312 int i, j, jmx, jj, jj0, k, k1, k2, l, count, numR[MAXV];
313 double fctr;
314 aNode *nd;
315 VV *vx, *vx1;
316 SS *sm;
317 Re *re;
318 Gem *gm;
319
320 /* some i/o */
321 Vnm_print(0,"Aprx_initNodes: initializing node structures..");
322
323 /* local pointers */
324 gm = thee->gm;
325
326 /* create the geometry map datastructures (empty) */
327 thee->numV = Gem_numVirtVV(gm);
328 thee->numE = Gem_numVirtEE(gm);
329 thee->numF = Gem_numVirtFF(gm);
330 thee->numS = Gem_numVirtSS(gm);
331 thee->I_V2E = Vmem_malloc( thee->vmem, 2*(thee->numE), sizeof(int) );
332 if (Gem_dim(gm) == 2) {
333 thee->I_V2S = Vmem_malloc( thee->vmem, 3*(thee->numS), sizeof(int) );
334 thee->I_E2S = Vmem_malloc( thee->vmem, 3*(thee->numS), sizeof(int) );
335 thee->I_V2F = VNULL;
336 thee->I_E2F = VNULL;
337 thee->I_F2S = VNULL;
338 } else if (Gem_dim(gm) == 3) {
339 thee->I_V2S = Vmem_malloc( thee->vmem, 4*(thee->numS), sizeof(int) );
340 thee->I_E2S = Vmem_malloc( thee->vmem, 6*(thee->numS), sizeof(int) );
341 thee->I_V2F = Vmem_malloc( thee->vmem, 3*(thee->numF), sizeof(int) );
342 thee->I_E2F = Vmem_malloc( thee->vmem, 3*(thee->numF), sizeof(int) );
343 thee->I_F2S = Vmem_malloc( thee->vmem, 4*(thee->numS), sizeof(int) );
344 } else { VASSERT(0); }
345
346 /* create the Bnode object (destroy any existing one) */
347 if (thee->B != VNULL) {
348 Bnode_dtor( &(thee->B) );
349 }
350 for (i=0; i<thee->numB; i++) {
351 re = thee->re[i];
352 numR[i] = Re_numVDF(re)*Gem_numVV(gm)
353 + Re_numEDF(re)*Gem_numVirtEE(gm)
354 + Re_numFDF(re)*Gem_numVirtVV(gm)
355 + Re_numSDF(re)*Gem_numSS(gm);
356 }
357 thee->B = Bnode_ctor(thee->vmem, thee->numB, numR);
358
359 /* Now initialize the Bnode */
360 count = 0;
361 for (i=0; i<thee->numB; i++) {
362 re = thee->re[i];
363 nd = Bnode_data(thee->B,i);
364 jmx = Re_numVDF(re)*Gem_numVV(gm);
365 for (j=0; j<jmx; j++) {
366 vx = Gem_VV(gm,j);
367 nd[j].myid = count;
368 nd[j].myvert = j; /* VV_id(vx); */
369 nd[j].type = 0;
370 nd[j].chart = VV_chart(vx);
371 nd[j].numx = Gem_dimII(gm);
372 nd[j].x[0] = VV_coord(vx,0);
373 nd[j].x[1] = VV_coord(vx,1);
374 nd[j].x[2] = VV_coord(vx,2);
375 nd[j].val = 0.0;
376 count++;
377 }
378
379 #ifdef VG_ELEMENT
380 /* Initialize edge degrees of freedom */
381 if (Re_numEDF(re)>0) {
382 jj0 = Re_numVDF(re)*Gem_numVirtVV(gm);
383 for (j=0; j<Gem_numSS(gm); j++) {
384 sm = Gem_SS(gm, j);
385 /* run through the edges */
386 for (k=0; k<Gem_dimEE(gm);k++) {
387 /* run through df on the edge */
388 for (l=0; l<Re_numEDF(re); l++) {
389 jj = l + Re_numEDF(re)*SS_edgeNumber(sm,k) + jj0;
390 if (nd[jj].numx != Gem_dimII(gm)) {
391 k1 = vmapEI[k][0];
392 k2 = vmapEI[k][1];
393 vx = SS_vertex (sm, k1);
394 vx1 = SS_vertex (sm, k2);
395 nd[jj].myid = count;
396 nd[jj].myvert = jj;
397 nd[jj].type = 0;
398 nd[jj].chart = 0;
399 nd[jj].numx = Gem_dimII(gm);
400 fctr = (double)(l+1)/(double)(Re_numEDF(re)+1);
401 nd[jj].x[0] = fctr * (VV_coord(vx,0)+VV_coord(vx1,0));
402 nd[jj].x[1] = fctr * (VV_coord(vx,1)+VV_coord(vx1,1));
403 nd[jj].x[2] = fctr * (VV_coord(vx,2)+VV_coord(vx1,2));
404 nd[jj].val = 0.0;
405 count++;
406 }
407 } /* l */
408 } /* k */
409 } /* j */
410 } /* if */
411 #endif
412
413 }
414
415 /* build new elements vectors for this mesh */
416 if (thee->wev != VNULL) {
417 Bvec_dtor( &(thee->wev) );
418 }
419 numR[0] = Gem_numSS(thee->gm);
420 thee->wev = Bvec_ctor( thee->vmem, "wev", 1, numR );
421
422 /* some i/o */
423 Vnm_print(0,"..done.\n");
424 }
425
426 /*
427 * ***************************************************************************
428 * Routine: Aprx_buildNodes
429 *
430 * Purpose: Determine the node types via a single fast traversal of
431 * the elements by looking at the element faces.
432 *
433 * Author: Michael Holst
434 * ***************************************************************************
435 */
Aprx_buildNodes(Aprx * thee)436 VPUBLIC void Aprx_buildNodes(Aprx *thee)
437 {
438 /* some variables */
439 int i, j, k, l, ef, es, ii, kk, kk0, numEF;
440 int vType, eType, fType, numBVD, vxid, count, chart;
441 double UUd[MAXV], xm[3];
442 SS *sm;
443 VV *vx;
444 Re *re;
445 Gem *gm;
446
447 /* initialization */
448 Aprx_initNodes(thee);
449 gm = thee->gm;
450
451 /* some i/o */
452 Vnm_print(0,"Aprx_buildNodes: determing node types from face types..");
453
454 /* dirichlet node count in each block */
455 for (ii=0; ii<thee->numB; ii++) {
456 numBVD = 0;
457 re = thee->re[ii];
458
459 /* fly through the elements */
460 for (i=0; i<Gem_numSS(gm); i++) {
461 sm = Gem_SS(gm,i);
462
463 /* fly through the faces of the element */
464 for (j=0; j<Gem_dimVV(gm); j++) {
465
466 /* here the user-defined boundary type remapping is used! */
467 fType = thee->pde->bmap[ii][ SS_faceType(sm,j) ];
468
469 /* if boundary face, apply type to vertices of the face */
470 if (VBOUNDARY( fType )) {
471
472 #ifdef VG_ELEMENT
473 /* type edge degrees of freedom */
474 if (Re_numEDF(re)>0) {
475
476 /* compute number of edges for a simplex 1 dimension lower */
477 numEF = (Re_dimV(re)-1)*(Re_dimV(re)-2)/2;
478 kk0 = Re_numVDF(re) * Gem_numVirtVV(gm);
479
480 /* run through the edges of the face */
481 for (ef=0; ef<numEF; ef++) {
482
483 /* get edge number wrt simplex */
484 es = vmapFE[j][ef];
485
486 /* run through df on the edge */
487 for (l=0; l<Re_numEDF(re); l++) {
488 k = l + Re_numEDF(re)*SS_edgeNumber(sm,es) + kk0;
489 eType = Bnode_data(thee->B,ii)[k].type;
490
491 /* if face is dirichlet, type is "hard" */
492 if (VDIRICHLET( fType )) {
493
494 /* is this a new dirichlet addition */
495 if ( ! VDIRICHLET( eType ) ) {
496 numBVD++;
497 }
498
499 /* type this edge */
500 Bnode_data(thee->B,ii)[k].type = fType;
501
502 /* if face is neumann, type is "soft" */
503 } else if (VNEUMANN( fType )) {
504
505 /* type this edge */
506 if ( ! VDIRICHLET( eType ) ) {
507 Bnode_data(thee->B,ii)[k].type = fType;
508 }
509
510 } else { VASSERT(0); }
511 } /* l */
512 } /* ef */
513 } /* if */
514 #endif
515
516 /* go through vertices of this face */
517 k=j;
518 for (kk=0; kk<Gem_dim(gm); kk++) {
519
520 /* get the vertex */
521 k=(k+1) % Gem_dimVV(gm);
522 vx = SS_vertex(sm,k);
523 vxid = VV_id(vx);
524 vType = Bnode_data(thee->B,ii)[vxid].type;
525
526 /* if face is dirichlet, type is "hard" */
527 if (VDIRICHLET( fType )) {
528
529 /* is this a new dirichlet addition */
530 if ( ! VDIRICHLET( vType ) ) {
531 numBVD++;
532 }
533
534 /* type this vertex */
535 Bnode_data(thee->B,ii)[vxid].type = fType;
536
537 /* if face is neumann, type is "soft" */
538 } else if (VNEUMANN( fType )) {
539
540 /* type this vertex */
541 if ( ! VDIRICHLET( vType ) ) {
542 Bnode_data(thee->B,ii)[vxid].type = fType;
543 }
544
545 } else { VASSERT(0); }
546 }
547 }
548 }
549 }
550
551 /* record node count, create space for dirichlet vertex list info */
552 if (thee->ibv[ii] != VNULL) {
553 Vmem_free(thee->vmem,
554 thee->numBV[ii],sizeof(int),(void**)&(thee->ibv[ii]));
555 }
556 if (thee->bv[ii] != VNULL) {
557 Vmem_free(thee->vmem,
558 thee->numBV[ii],sizeof(double),(void**)&(thee->bv[ii]));
559 }
560 thee->numBV[ii] = numBVD;
561 if (thee->numBV[ii] > 0) {
562 thee->ibv[ii] = Vmem_malloc(thee->vmem,
563 thee->numBV[ii],sizeof(int));
564 thee->bv[ii] = Vmem_malloc(thee->vmem,
565 thee->numBV[ii],sizeof(double));
566 } else {
567 thee->ibv[ii] = VNULL;
568 thee->bv[ii] = VNULL;
569 }
570 }
571
572 /* initalize dirichlet boundary vertex list */
573 for (k=0; k<thee->numB; k++) {
574 count = 0;
575 for (i=0; i<Bnode_numR(thee->B,k); i++) {
576 for (j=0; j<Bnode_data(thee->B,k)[i].numx; j++) {
577 xm[j] = Bnode_data(thee->B,k)[i].x[j];
578 }
579 chart = Bnode_data(thee->B,k)[i].chart;
580 vType = Bnode_data(thee->B,k)[i].type;
581 thee->pde->u_D(thee->pde, vType, chart, xm, UUd);
582 if (VDIRICHLET( Bnode_data(thee->B,k)[i].type )) {
583 Bnode_data(thee->B,k)[i].val = UUd[k];
584 thee->ibv[k][count] = i;
585 thee->bv[k][count] = UUd[k];
586 count++;
587 }
588 }
589 Vnm_print(0,"count = %d, thee->numBV[0] = %d\n", count, thee->numBV[k]);
590 VASSERT( count == thee->numBV[k] );
591 }
592
593 /* some i/o */
594 Vnm_print(0,"..done.\n");
595 }
596
597 /*
598 * ***************************************************************************
599 * Routine: Aprx_evalTrace
600 *
601 * Purpose: Evaluate the trace information.
602 *
603 * Author: Michael Holst
604 * ***************************************************************************
605 */
Aprx_evalTrace(Aprx * thee,Bvec * Wud,Bvec * Wui,Bvec * Wut)606 VPUBLIC void Aprx_evalTrace(Aprx *thee, Bvec *Wud, Bvec *Wui, Bvec *Wut)
607 {
608 /* some variables */
609 int i, j, k, vType, count, chart;
610 double UUd[MAXV], UUt[MAXV], dUUt[MAXV][3], xm[3];
611
612 /* some i/o */
613 Vnm_print(0,"Aprx_evalTrace: evaluating <Ud,Ui,Ut>..");
614
615 /* initalize pde solution and dirichlet boundary vertex list */
616 for (k=0; k<thee->numB; k++) {
617 count = 0;
618 for (i=0; i<Bnode_numR(thee->B,k); i++) {
619 for (j=0; j<Bnode_data(thee->B,k)[i].numx; j++) {
620 xm[j] = Bnode_data(thee->B,k)[i].x[j];
621 }
622 chart = Bnode_data(thee->B,k)[i].chart;
623 vType = Bnode_data(thee->B,k)[i].type;
624 thee->pde->u_D(thee->pde, vType, chart, xm, UUd);
625 thee->pde->u_T(thee->pde, vType, chart, xm, UUt, dUUt);
626 if (VDIRICHLET( Bnode_data(thee->B,k)[i].type )) {
627 Bvec_setB( Wud, k, i, UUd[k] );
628 Bvec_setB( Wui, k, i, 0. );
629 Bvec_setB( Wut, k, i, 0. );
630 count++;
631 } else {
632 Bvec_setB( Wud, k, i, 0. );
633 Bvec_setB( Wui, k, i, UUd[k] );
634 Bvec_setB( Wut, k, i, UUt[k] );
635 }
636 }
637 VASSERT( count == thee->numBV[k] );
638 }
639
640 /* some i/o */
641 Vnm_print(0,"..done.\n");
642 }
643
644 /*
645 * ***************************************************************************
646 * Routine: Aprx_refine
647 *
648 * Purpose: Refine the mesh.
649 *
650 * Author: Michael Holst
651 * ***************************************************************************
652 */
Aprx_refine(Aprx * thee,int rkey,int bkey,int pkey)653 VPUBLIC int Aprx_refine(Aprx *thee, int rkey, int bkey, int pkey)
654 {
655 int numRc, numR, refined;
656
657 /* number of vertices before refinement */
658 numRc = Gem_numVV(thee->gm);
659
660 /* call refinement */
661 refined = Gem_refine(thee->gm, rkey, bkey, pkey);
662
663 /* number of vertices after refinement */
664 numR = Gem_numVV(thee->gm);
665
666 /* use lazy edge creation to build a linear prolongation operator */
667 Aprx_buildProlong(thee, numRc, numR, pkey);
668
669 /* destroy lazy edges that we used to build prolongation operator */
670 Gem_destroyEdges(thee->gm);
671
672 /* count v/e/f/s and check the mesh */
673 Gem_countChk(thee->gm);
674
675 /* return */
676 return refined;
677 }
678
679 /*
680 * ***************************************************************************
681 * Routine: Aprx_unRefine
682 *
683 * Purpose: Unrefine the mesh.
684 *
685 * Author: Michael Holst
686 * ***************************************************************************
687 */
Aprx_unRefine(Aprx * thee,int rkey,int pkey)688 VPUBLIC int Aprx_unRefine(Aprx *thee, int rkey, int pkey)
689 {
690 int refined;
691
692 /* call unrefinement */
693 refined = Gem_unRefine(thee->gm, rkey, pkey);
694
695 /* return */
696 return refined;
697 }
698
699 /*
700 * ***************************************************************************
701 * Routine: Aprx_buildProlong
702 *
703 * Purpose: Create storage for and build prolongation operator; the operator
704 * prolongates a vector FROM a coarser level TO this level.
705 *
706 * Notes: This routine assumes the prolongation operator, P = [ I ],
707 * where I is the identity and R is the tail of P. [ R ]
708 * For "single-refined" edges, the rows of R are
709 * zero, except for 1/2 in the columns corresponding to the
710 * parent vertices. For a multiply refined edge, we just
711 * take a half-half linear combination of the parent rows.
712 * We first build R as an RLN, then copy it into P in ROW format.
713 *
714 * Author: Stephen Bond and Michael Holst
715 * ***************************************************************************
716 */
Aprx_buildProlong(Aprx * thee,int numRc,int numR,int pkey)717 VPUBLIC void Aprx_buildProlong(Aprx *thee, int numRc, int numR, int pkey)
718 {
719 int k, i, j0, j1, count, queuecount, multRefCount;
720 int pnumR, pnumC, pnumP, ibase, *IJP;
721 EE *eg;
722 VV *v01, *v0, *v1;
723 Mat *Pmat, *Rmat;
724 LinkA *mt, *mt0, *mt1;
725
726 Vnm_print(0,"Aprx_buildProlong: creating prolongation operator..");
727
728 /* Determine dimensions of new prolongation operator */
729 pnumR = numR;
730 pnumC = numRc;
731 ibase = (pkey) ? pnumC : 0;
732
733 /* Construct a matrix for the tail of the new prolongation operator */
734 Rmat = Mat_ctor( VNULL, "P", (pnumR-pnumC), pnumC );
735
736 /* Initialize Rmat with the RLN format */
737 Mat_initStructureLN( Rmat, RLN_FORMAT, ISNOT_SYM );
738
739 /* STEP 1: Loop over split edges, filling in the entries of Rmat */
740 count = 0;
741 queuecount = 0;
742 multRefCount = 0;
743 for (k=0; k<Gem_numEE(thee->gm); k++) {
744 eg = Gem_EE(thee->gm,k);
745
746 /* Get the new vertex */
747 v01 = EE_midPoint(eg);
748 VASSERT( v01 != VNULL );
749 i = VV_id(v01) - pnumC;
750 VASSERT( i >= 0 );
751
752 /* Get the parent vertices */
753 v0 = EE_vertex(eg,0);
754 j0 = VV_id(v0);
755 v1 = EE_vertex(eg,1);
756 j1 = VV_id(v1);
757
758 /* Get the head of the current row */
759 mt = (LinkA*)Vset_access(Rmat->lnkU,i);
760
761 /* Have we already completed this row? */
762 if( mt->idx < 0 ) {
763
764 /* Case 1: Vertex v0 is old and Vertex v1 is old */
765 if( (j0 < pnumC) && (j1 < pnumC) ) {
766 mContrib( Rmat->lnkU, 1, &count, i, j0, 0.5 );
767 mContrib( Rmat->lnkU, 1, &count, i, j1, 0.5 );
768 }
769
770 /* Case 2: Vertex v0 is new and Vertex v1 is old */
771 if( (j0 >= pnumC) && (j1 < pnumC) ) {
772 mt0 = (LinkA*)Vset_access(Rmat->lnkU,j0-pnumC);
773 if( mt0->idx < 0 ) {
774 /* Parent row unresolved */
775 queuecount++;
776 } else {
777 /* Parent row resolved */
778 multRefCount++;
779 mContrib( Rmat->lnkU, 1, &count, i, j1, 0.5 );
780 for (mt=mt0; mt!=VNULL; mt=mt->next) {
781 mContrib( Rmat->lnkU, 1, &count,
782 i, mt->idx, 0.5*(mt->val) );
783 }
784 }
785 }
786
787 /* Case 3: Vertex v0 is old and Vertex v1 is new */
788 if( (j0 < pnumC) && (j1 >= pnumC) ) {
789 mt1 = (LinkA*)Vset_access(Rmat->lnkU,j1-pnumC);
790 if( mt1->idx < 0 ) {
791 /* Parent row unresolved */
792 queuecount++;
793 } else {
794 /* Parent row resolved */
795 multRefCount++;
796 mContrib( Rmat->lnkU, 1, &count, i, j0, 0.5 );
797 for (mt=mt1; mt!=VNULL; mt=mt->next) {
798 mContrib( Rmat->lnkU, 1, &count,
799 i, mt->idx, 0.5*(mt->val) );
800 }
801 }
802 }
803
804 /* Case 4: Vertex v0 is new and Vertex v1 is new */
805 if( (j0 >= pnumC) && (j1 >= pnumC) ) {
806 mt0 = (LinkA*)Vset_access(Rmat->lnkU,j0-pnumC);
807 mt1 = (LinkA*)Vset_access(Rmat->lnkU,j1-pnumC);
808 if( mt0->idx < 0 || mt1->idx < 0 ) {
809 /* Parent row(s) unresolved */
810 queuecount += 2;
811 } else {
812 /* Parent row(s) resolved */
813 multRefCount++;
814 for (mt=mt0; mt!=VNULL; mt=mt->next) {
815 mContrib( Rmat->lnkU, 1, &count,
816 i, mt->idx, 0.5*(mt->val) );
817 }
818 multRefCount++;
819 for (mt=mt1; mt!=VNULL; mt=mt->next) {
820 mContrib( Rmat->lnkU, 1, &count,
821 i, mt->idx, 0.5*(mt->val) );
822 }
823 }
824 }
825 }
826 }
827 if (multRefCount > 0) {
828 Vnm_print(0,"[RESOLVED MULTIPLY REFINED EDGES = <%d>]",multRefCount);
829 }
830 if (queuecount > 0) {
831 /* ERROR: This should never happen if edges are processed in order! */
832 Vnm_print(0,"[UNRESOLVED MULTIPLY REFINED EDGES = <%d>!]",queuecount);
833 /* We could add a loop until queue is empty, but is not necessary? */
834 VASSERT(0);
835 }
836
837 /* STEP 2: Construct Pmat from Rmat */
838 Pmat = Mat_ctor( VNULL, "P", (pnumR-ibase), pnumC );
839
840 /* create actual space for new prolongation operator */
841 pnumP = (pnumC - ibase) + count;
842
843 IJP = Vmem_malloc( Pmat->vmem, (pnumR-ibase)+1+pnumP, sizeof(int) );
844 Mat_initStructure( Pmat, ROW_FORMAT, ISNOT_SYM, pnumP, IJP, VNULL );
845 Mat_setState( Pmat, ASSEMBLED_STATE );
846 if (pkey == 1) {
847 Mat_setImpl( Pmat, IS_IMPL );
848 }
849
850 /* STEP 2A: Create the identity for the shared coarse and fine points */
851 count = 0;
852 Pmat->IA[count] = 0;
853 for (i=0; i<(pnumC-ibase); i++) {
854 Pmat->A[count] = 1.;
855 Pmat->JA[count] = i;
856 count++;
857 Pmat->IA[i+1] = count;
858 }
859
860 /* STEP 2B: Copy Rmat (RLN) into P (ROW) */
861 for (i=0; i<Rmat->numR; i++) {
862 for (mt=(LinkA*)Vset_access(Rmat->lnkU,i); mt!=VNULL; mt=mt->next) {
863 if( mt->idx >= 0 ) {
864 Pmat->A[count] = mt->val;
865 Pmat->JA[count] = mt->idx;
866 count++;
867 }
868 }
869 Pmat->IA[(pnumC-ibase)+i+1] = count;
870 }
871
872 /* Clean-up: Destroy temporary Rmat */
873 Mat_dtor( &Rmat );
874
875 Vnm_print(0,"..done.\n");
876
877 /* STEP 3: Convert to block version (also updates Bnode) */
878 Aprx_buildProlongBmat(thee, Pmat);
879 Mat_dtor( &Pmat );
880 }
881
882 /*
883 * ***************************************************************************
884 * Routine: Aprx_buildProlongBmat
885 *
886 * Purpose: Create block version of prolongation operator; also
887 * updates the Bnode object.
888 *
889 * Author: Michael Holst
890 * ***************************************************************************
891 */
Aprx_buildProlongBmat(Aprx * thee,Mat * p)892 VPUBLIC void Aprx_buildProlongBmat(Aprx *thee, Mat *p)
893 {
894 Bmat *bp;
895 int i, j, numR[MAXV], numC[MAXV];
896 MATmirror mirror[MAXV][MAXV];
897 int numBR[MAXV], *BR[MAXV];
898 int numBC[MAXV], *BC[MAXV];
899
900 /* prolongation matrix structure check */
901 VASSERT(p != VNULL);
902 VASSERT( Mat_format(p) == ROW_FORMAT );
903
904 /* setup mirror information */
905 for (i=0; i<thee->numB; i++) {
906 numR[i] = Mat_numR(p);
907 numC[i] = Mat_numC(p);
908 for (j=0; j<thee->numB; j++) {
909 mirror[i][j] = ISNOT_MIRROR;
910 }
911 }
912
913 /* copy old boundary node arrays while we still have them */
914 for (i=0; i<thee->numB; i++) {
915 numBC[i] = thee->numBV[i];
916 if (numBC[i] > 0) {
917 BC[i] = Vmem_malloc(thee->vmem,numBC[i],sizeof(int));
918 }
919 for (j=0; j<numBC[i]; j++) {
920 BC[i][j] = thee->ibv[i][j];
921 }
922 }
923
924 /* now build new nodes */
925 Aprx_buildNodes(thee);
926
927 /* now copy new boundary node arrays */
928 for (i=0; i<thee->numB; i++) {
929 numBR[i] = thee->numBV[i];
930 if (numBR[i] > 0) {
931 BR[i] = Vmem_malloc(thee->vmem,numBR[i],sizeof(int));
932 }
933 for (j=0; j<numBR[i]; j++) {
934 BR[i][j] = thee->ibv[i][j];
935 }
936 }
937
938 /* build new prolongation; keep pointer to old "coarse" prolongation */
939 bp = Bmat_ctor( thee->vmem, "P", thee->numB, numR, numC, mirror );
940 for (i=0; i<thee->numB; i++) {
941
942 /* copy the nonzero structure */
943 Mat_copyStructure( bp->AD[i][i], p );
944
945 /* copy the nonzeros themselves and mark as assembled */
946 for (j=0; j<p->numA; j++) {
947 bp->AD[i][i]->A[j] = p->A[j];
948 }
949 Mat_setState( bp->AD[i][i], ASSEMBLED_STATE );
950
951 /* apply any dirichlet node conditions */
952 Mat_buildBRC(bp->AD[i][i], numBR[i], numBC[i], BR[i], BC[i]);
953 Mat_zeroBRC(bp->AD[i][i]);
954 }
955
956 /* release local boundary node arrays */
957 for (i=0; i<thee->numB; i++) {
958 if (numBR[i] > 0) {
959 Vmem_free(thee->vmem,numBR[i],sizeof(int),(void**)&BR[i]);
960 }
961 if (numBC[i] > 0) {
962 Vmem_free(thee->vmem,numBC[i],sizeof(int),(void**)&BC[i]);
963 }
964 }
965
966 /* setup the prolongation chain */
967 bp->coarse = thee->P;
968 thee->P = bp;
969 }
970
971 /*
972 * ***************************************************************************
973 * Routine: bContrib
974 *
975 * Purpose: Set or add an entry to a linked matrix entry list.
976 *
977 * Author: Michael Holst
978 * ***************************************************************************
979 */
bContrib(Vset * mtpool,int * count,int * countyr,int i,int j)980 VPRIVATE void bContrib(Vset *mtpool, int *count, int *countyr, int i, int j)
981 {
982 int done;
983 Link *curr, *mnew;
984
985 mnew=VNULL;
986 curr=(Link*)Vset_access(mtpool,i);
987 VASSERT( curr != VNULL );
988
989 /* we have to look through the row(col) */
990 done=0;
991 while (!done) {
992
993 /* if first guy in row(col) is "blank", fill him with the info */
994 if (curr->idx == -1) {
995 (*countyr)++;
996 if (i<j) (*count)++;
997 curr->idx = j;
998 curr->next = VNULL;
999 done = 1;
1000
1001 /* we found the position; already here so nothing to do */
1002 } else if (curr->idx == j) {
1003 done = 1;
1004
1005 /* not in list; insert new struct AFTER him, then swap me and him */
1006 } else if (curr->idx > j) {
1007 (*countyr)++;
1008 if (i<j) (*count)++;
1009 mnew = (Link*)Vset_create(mtpool);
1010 mnew->idx = curr->idx;
1011 mnew->next = curr->next;
1012 curr->idx = j;
1013 curr->next = mnew;
1014 done = 1;
1015
1016 /* not found; no one left so create a new guy and put me there */
1017 } else if (curr->next == VNULL) {
1018 (*countyr)++;
1019 if (i<j) (*count)++;
1020 mnew = (Link*)Vset_create(mtpool);
1021 mnew->idx = j;
1022 mnew->next = VNULL;
1023 curr->next = mnew;
1024 done = 1;
1025
1026 /* not found yet; still hope */
1027 } else {
1028 curr=curr->next;
1029 }
1030 }
1031 }
1032
1033 /*
1034 * ***************************************************************************
1035 * Routine: Aprx_nodeCount
1036 *
1037 * Purpose: Count up the total number of basis functions (nodes).
1038 * This is basically an a priori count of the number of rows in
1039 * the stiffness matrix.
1040 *
1041 * Notes: This routine can handle COMPLETELY GENERAL ELEMENTS; see the
1042 * comments in Aprx_interact().
1043 *
1044 * Author: Michael Holst
1045 * ***************************************************************************
1046 */
Aprx_nodeCount(Aprx * thee,Re * re)1047 VPUBLIC int Aprx_nodeCount(Aprx *thee, Re *re)
1048 {
1049 int nodeCount;
1050
1051 /* make sure we have Vertex/Edge/Face/Simplex counts if needed */
1052 if (Re_numVDF(re) > 0) VASSERT( Gem_numVirtVV(thee->gm) > 0 );
1053 if (Re_numEDF(re) > 0) VASSERT( Gem_numVirtEE(thee->gm) > 0 );
1054 if (Re_numFDF(re) > 0) VASSERT( Gem_numVirtFF(thee->gm) > 0 );
1055 if (Re_numSDF(re) > 0) VASSERT( Gem_numVirtSS(thee->gm) > 0 );
1056
1057 /* calculate the node count */
1058 nodeCount = Gem_numVirtVV(thee->gm) * Re_numVDF(re)
1059 + Gem_numVirtEE(thee->gm) * Re_numEDF(re)
1060 + Gem_numVirtFF(thee->gm) * Re_numFDF(re)
1061 + Gem_numVirtSS(thee->gm) * Re_numSDF(re);
1062
1063 /* return */
1064 return nodeCount;
1065 }
1066
1067 /*
1068 * ***************************************************************************
1069 * Routine: Aprx_interact
1070 *
1071 * Purpose: Count the number of basis (linear/quadratic/etc) functions which
1072 * interact with a given basis function. This is basically an
1073 * a priori count of the number of nonzeros per row in the
1074 * stiffness matrix, and the number of nozeros per row above the
1075 * diagonal (and thus the number of nonzeros per column below the
1076 * diagonal if we have symmetric nonzero structure). By doing this
1077 * count, we can one-time malloc storage for matrices before
1078 * assembly, and avoid extra copies as well as avoid the use
1079 * of linked list structures for representing matrices with
1080 * a priori unknown nonzero structure.
1081 *
1082 * Notes: This routine can handle COMPLETELY GENERAL ELEMENTS; all it needs
1083 * to do is simply compute all interactions of all nodes in the
1084 * element. The nodes may be any combination of vertex-based,
1085 * edge-based, face-based (3D only), or simplex-based (interior)
1086 * degrees of freedom.
1087 *
1088 * This calculation is made possible because all of the node numbers
1089 * for all of the vertex/edge/face/simplex-nodes in each simplex are
1090 * available in O(1) time for each element, using the geometry
1091 * manager routine Gem_simplexInfo().
1092 *
1093 * Author: Michael Holst (updated by Ryan Szypowski)
1094 * ***************************************************************************
1095 */
Aprx_interact(Aprx * thee,Re * re,Re * reT,int * numR,int * numO,int * numOYR,int ** IJA,int ** IJAYR)1096 VPUBLIC void Aprx_interact(Aprx *thee, Re *re, Re *reT,
1097 int *numR, int *numO, int *numOYR, int **IJA, int **IJAYR)
1098 {
1099 int smid, i, j, p, q, ii, jj, count, countyr, size;
1100 int *IA,*JA,*IAYR,*JAYR;
1101 TT tt;
1102 Link *mt;
1103 SS *sm;
1104 Gem *gm;
1105
1106 /* go through all nodes and setup the node interaction structure */
1107 Vnm_tstart(40, "node interactions");
1108 Vnm_print(0,"Aprx_interact: counting node interactions..\n");
1109
1110 /* some local pointers */
1111 gm = thee->gm;
1112
1113 /* clear the dynamic array */
1114 Vset_reset( thee->lnkG );
1115
1116 /* calculate some node information */
1117 *numR = Aprx_nodeCount(thee, re);
1118
1119 /* create an empty entry to start each row of global matrix */
1120 for (i=0; i<*numR; i++) {
1121 mt=(Link*)Vset_create(thee->lnkG);
1122 mt->idx = -1;
1123 mt->next = VNULL;
1124 }
1125
1126 /* go through simplices and count the off-diagonal interactions */
1127 count = 0;
1128 countyr = 0;
1129
1130 /* offsets for edge/face/simplex DF global indices */
1131 re->offIS[0] = 0;
1132 re->offIS[1] = Re_numVDF(re) * Gem_numVirtVV(gm);
1133 re->offIS[2] = Re_numEDF(re) * Gem_numVirtEE(gm) + re->offIS[1];
1134 re->offIS[3] = Re_numFDF(re) * Gem_numVirtFF(gm) + re->offIS[2];
1135 #if 0
1136 re->offIS[4] = Re_numSDF(re) * Gem_numVirtSS(gm) + re->offIS[3];
1137 #endif
1138
1139 reT->offIS[0] = 0;
1140 reT->offIS[1] = Re_numVDF(reT) * Gem_numVirtVV(gm);
1141 reT->offIS[2] = Re_numEDF(reT) * Gem_numVirtEE(gm) + reT->offIS[1];
1142 reT->offIS[3] = Re_numFDF(reT) * Gem_numVirtFF(gm) + reT->offIS[2];
1143 #if 0
1144 reT->offIS[4] = Re_numSDF(reT) * Gem_numVirtSS(gm) + reT->offIS[3];
1145 #endif
1146
1147 for (smid=0; smid<Gem_numSS(gm); smid++) {
1148 sm = Gem_SS(gm,smid);
1149 Gem_simplexInfo(gm, sm, &tt);
1150
1151 /* calculate interactions: vertex-* */
1152 for (p=0; p<Re_numVDF(re); p++) {
1153 for (i=0; i<tt.dimV; i++) {
1154 ii = Aprx_nodeIndex (thee, re, tt.vid[i], 0, p);
1155
1156 /* calculate interactions: vertex-vertex */
1157 for (q=0; q<Re_numVDF(reT); q++) {
1158 for (j=0; j<tt.dimV; j++) {
1159 jj = Aprx_nodeIndex (thee, reT, tt.vid[j], 0, q);
1160 bContrib(thee->lnkG,&count,&countyr,ii,jj);
1161 }
1162 }
1163
1164 /* calculate interactions: vertex-edge */
1165 for (q=0; q<Re_numEDF(reT); q++) {
1166 for (j=0; j<tt.dimE; j++) {
1167 jj = Aprx_nodeIndex (thee, reT, tt.eid[j], 1, q);
1168 bContrib(thee->lnkG,&count,&countyr,ii,jj);
1169 }
1170 }
1171
1172 /* calculate interactions: vertex-face */
1173 for (q=0; q<Re_numFDF(reT); q++) {
1174 for (j=0; j<tt.dimF; j++) {
1175 jj = Aprx_nodeIndex (thee, reT, tt.fid[j], 2, q);
1176 bContrib(thee->lnkG,&count,&countyr,ii,jj);
1177 }
1178 }
1179
1180 /* calculate interactions: vertex-simplex */
1181 for (q=0; q<Re_numSDF(reT); q++) {
1182 jj = Aprx_nodeIndex (thee, reT, tt.sid, 3, q);
1183 bContrib(thee->lnkG,&count,&countyr,ii,jj);
1184 }
1185 }
1186 }
1187
1188 /* calculate interactions: edge-* */
1189 for (p=0; p<Re_numEDF(re); p++) {
1190 for (i=0; i<tt.dimE; i++) {
1191 ii = Aprx_nodeIndex (thee, re, tt.eid[i], 1, p);
1192
1193 /* calculate interactions: edge-vertex */
1194 for (q=0; q<Re_numVDF(reT); q++) {
1195 for (j=0; j<tt.dimV; j++) {
1196 jj = Aprx_nodeIndex (thee, reT, tt.vid[j], 0, q);
1197 bContrib(thee->lnkG,&count,&countyr,ii,jj);
1198 }
1199 }
1200
1201 /* calculate interactions: edge-edge */
1202 for (q=0; q<Re_numEDF(reT); q++) {
1203 for (j=0; j<tt.dimE; j++) {
1204 jj = Aprx_nodeIndex (thee, reT, tt.eid[j], 1, q);
1205 bContrib(thee->lnkG,&count,&countyr,ii,jj);
1206 }
1207 }
1208
1209 /* calculate interactions: edge-face */
1210 for (q=0; q<Re_numFDF(reT); q++) {
1211 for (j=0; j<tt.dimF; j++) {
1212 jj = Aprx_nodeIndex (thee, reT, tt.fid[j], 2, q);
1213 bContrib(thee->lnkG,&count,&countyr,ii,jj);
1214 }
1215 }
1216
1217 /* calculate interactions: edge-simplex */
1218 for (q=0; q<Re_numSDF(reT); q++) {
1219 jj = Aprx_nodeIndex (thee, reT, tt.sid, 3, q);
1220 bContrib(thee->lnkG,&count,&countyr,ii,jj);
1221 }
1222 }
1223 }
1224
1225 /* calculate interactions: face-* */
1226 for (p=0; p<Re_numFDF(re); p++) {
1227 for (i=0; i<tt.dimF; i++) {
1228 ii = Aprx_nodeIndex (thee, re, tt.fid[i], 2, p);
1229
1230 /* calculate interactions: face-vertex */
1231 for (q=0; q<Re_numVDF(reT); q++) {
1232 for (j=0; j<tt.dimV; j++) {
1233 jj = Aprx_nodeIndex (thee, reT, tt.vid[j], 0, q);
1234 bContrib(thee->lnkG,&count,&countyr,ii,jj);
1235 }
1236 }
1237
1238 /* calculate interactions: face-edge */
1239 for (q=0; q<Re_numEDF(reT); q++) {
1240 for (j=0; j<tt.dimE; j++) {
1241 jj = Aprx_nodeIndex (thee, reT, tt.eid[j], 1, q);
1242 bContrib(thee->lnkG,&count,&countyr,ii,jj);
1243 }
1244 }
1245
1246 /* calculate interactions: face-face */
1247 for (q=0; q<Re_numFDF(reT); q++) {
1248 for (j=0; j<tt.dimF; j++) {
1249 jj = Aprx_nodeIndex (thee, reT, tt.fid[j], 2, q);
1250 bContrib(thee->lnkG,&count,&countyr,ii,jj);
1251 }
1252 }
1253
1254 /* calculate interactions: face-simplex */
1255 for (q=0; q<Re_numSDF(reT); q++) {
1256 jj = Aprx_nodeIndex (thee, reT, tt.sid, 3, q);
1257 bContrib(thee->lnkG,&count,&countyr,ii,jj);
1258 }
1259 }
1260 }
1261
1262 /* calculate interactions: simplex-* */
1263 for (p=0; p<Re_numSDF(re); p++) {
1264 ii = Aprx_nodeIndex (thee, re, tt.sid, 3, p);
1265
1266 /* calculate interactions: simplex-vertex */
1267 for (q=0; q<Re_numVDF(reT); q++) {
1268 for (j=0; j<tt.dimV; j++) {
1269 jj = Aprx_nodeIndex (thee, reT, tt.vid[j], 0, q);
1270 bContrib(thee->lnkG,&count,&countyr,ii,jj);
1271 }
1272 }
1273
1274 /* calculate interactions: simplex-edge */
1275 for (q=0; q<Re_numEDF(reT); q++) {
1276 for (j=0; j<tt.dimE; j++) {
1277 jj = Aprx_nodeIndex (thee, reT, tt.eid[j], 1, q);
1278 bContrib(thee->lnkG,&count,&countyr,ii,jj);
1279 }
1280 }
1281
1282 /* calculate interactions: simplex-face */
1283 for (q=0; q<Re_numFDF(reT); q++) {
1284 for (j=0; j<tt.dimF; j++) {
1285 jj = Aprx_nodeIndex (thee, reT, tt.fid[j], 2, q);
1286 bContrib(thee->lnkG,&count,&countyr,ii,jj);
1287 }
1288 }
1289
1290 /* calculate interactions: simplex-simplex */
1291 for (q=0; q<Re_numSDF(reT); q++) {
1292 jj = Aprx_nodeIndex (thee, reT, tt.sid, 3, q);
1293 bContrib(thee->lnkG,&count,&countyr,ii,jj);
1294 }
1295 }
1296 }
1297 *numO = count;
1298 *numOYR = countyr;
1299 Vnm_print(0,"Aprx_interact: numN=<%d>, numO=<%d> numOYR=<%d>\n",
1300 *numR, *numO, *numOYR);
1301 VASSERT( *numOYR == (*numR+2*(*numO)) );
1302
1303 /* create the row-start and column-index arrays for YSMP-bank */
1304 size = *numR + 1 + *numO;
1305 *IJA = Vmem_malloc( thee->vmem, size, sizeof(int) );
1306 IA = *IJA;
1307 JA = *IJA + *numR + 1;
1308
1309 /* create the row-start and column-index arrays for YSMP-row(col) */
1310 size = *numR + 1 + *numOYR;
1311 *IJAYR = Vmem_malloc( thee->vmem, size, sizeof(int) );
1312 IAYR = *IJAYR;
1313 JAYR = *IJAYR + *numR + 1;
1314
1315 /* now create the matrix from the temporary structure */
1316 count = 0;
1317 countyr = 0;
1318 IA[0] = 0;
1319 IAYR[0] = 0;
1320 for (i=0; i<*numR; i++) {
1321 for (mt=(Link*)Vset_access(thee->lnkG,i); mt!=VNULL; mt=mt->next) {
1322 if (mt->idx >= 0) {
1323 JAYR[countyr] = mt->idx;
1324 countyr++;
1325 if (mt->idx > i) {
1326 JA[count] = mt->idx;
1327 count++;
1328 }
1329 }
1330 }
1331 IA[i+1] = count;
1332 IAYR[i+1] = countyr;
1333 }
1334 VASSERT( count == *numO );
1335 VASSERT( countyr == *numOYR );
1336
1337 /* clear the dynamic array */
1338 Vset_reset( thee->lnkG );
1339
1340 /* finish up */
1341 Vnm_print(0,"Aprx_interact: Done.\n");
1342 Vnm_tstop(40, "node interactions");
1343 }
1344
1345 /*
1346 * ***************************************************************************
1347 * Routine: Aprx_interactBlock
1348 *
1349 * Purpose: Determine all basis function interations.
1350 *
1351 * Author: Michael Holst (modified by Ryan Szypowski)
1352 * ***************************************************************************
1353 */
Aprx_interactBlock(Aprx * thee,Re * re[MAXV],Re * reT[MAXV],MATsym sym[MAXV][MAXV],MATmirror mirror[MAXV][MAXV],MATformat frmt[MAXV][MAXV],int numR[MAXV],int numO[MAXV][MAXV],int * IJA[MAXV][MAXV])1354 VPUBLIC void Aprx_interactBlock(Aprx *thee, Re *re[MAXV], Re *reT[MAXV],
1355 MATsym sym[MAXV][MAXV], MATmirror mirror[MAXV][MAXV], MATformat frmt[MAXV][MAXV],
1356 int numR[MAXV], int numO[MAXV][MAXV], int *IJA[MAXV][MAXV])
1357 {
1358 int i, j, k;
1359 int nr, no, noyr, *ija, *ijayr, *itmp;
1360
1361 for (i=0; i<thee->numB; i++) {
1362 for (j=0; j<thee->numB; j++) {
1363
1364 /* store the symmetry key for this block */
1365 switch (thee->pde->sym[i][j]) {
1366 case 0:
1367 sym[i][j] = ISNOT_SYM;
1368 mirror[i][j] = ISNOT_MIRROR;
1369 break;
1370 case 1:
1371 sym[i][j] = IS_SYM;
1372 mirror[i][j] = ISNOT_MIRROR;
1373 break;
1374 case 2:
1375 sym[i][j] = ISNOT_SYM; /* ignored when <mirror=IS_MIRROR> */
1376 mirror[i][j] = IS_MIRROR;
1377 break;
1378 default:
1379 VASSERT(0);
1380 break;
1381 }
1382
1383 /* upper-triangular blocks built first */
1384 if (j >= i) {
1385
1386 /* compute the interactions */
1387 Aprx_interact( thee, re[i], reT[j],
1388 &nr, &no, &noyr, &ija, &ijayr );
1389 numR[i] = nr;
1390
1391 /* handle diag/off-diag differences */
1392 if (j==i) {
1393
1394 /* diagonal blocks: bank-YSMP */
1395 frmt[i][j] = DRC_FORMAT;
1396 numO[i][j] = no;
1397 IJA[i][j] = ija;
1398 Vmem_free( thee->vmem, nr+1+noyr, sizeof(int),
1399 (void**)&ijayr );
1400
1401 } else if (j>i) {
1402
1403 /* upper-triangle blocks: row-YSMP */
1404 frmt[i][j] = ROW_FORMAT;
1405 numO[i][j] = noyr;
1406 IJA[i][j] = ijayr;
1407 Vmem_free( thee->vmem, nr+1+no, sizeof(int),
1408 (void**)&ija );
1409
1410 }
1411
1412 } else { /* (j<i) */
1413
1414 if ( !mirror[i][j] ) {
1415
1416 /* build lower-triangle block from upper triangle block */
1417 nr = numR[j];
1418 noyr = numO[j][i];
1419 ijayr = Vmem_malloc( thee->vmem, nr+1+noyr, sizeof(int) );
1420 itmp = IJA[j][i];
1421 for (k=0; k<(nr+1+noyr); k++) {
1422 ijayr[k] = itmp[k];
1423 }
1424
1425 /* lower-triangle blocks: col-YSMP */
1426 frmt[i][j] = COL_FORMAT;
1427 numO[i][j] = noyr;
1428 IJA[i][j] = ijayr;
1429
1430 }
1431 }
1432 }
1433 }
1434
1435 }
1436
1437 /*
1438 * ***************************************************************************
1439 * Routine: Aprx_nodeIndex
1440 *
1441 * Purpose: Compute global index of node of dimension dim in the matrix
1442 *
1443 * Author: Michael Holst, Oleg Korobkin (updated by Ryan Szypowski)
1444 * ***************************************************************************
1445 */
Aprx_nodeIndex(Aprx * thee,Re * re,int i,int dim,int q)1446 VPUBLIC int Aprx_nodeIndex (Aprx *thee, Re *re, int i, int dim, int q)
1447 {
1448 int ir;
1449
1450 ir = re->offIS[dim] + i;
1451 switch (dim) {
1452 case 0: ir += q * Gem_numVirtVV(thee->gm); break;
1453 case 1: ir += q * Gem_numVirtEE(thee->gm); break;
1454 case 2: ir += q * Gem_numVirtFF(thee->gm); break;
1455 case 3: ir += q * Gem_numVirtSS(thee->gm); break;
1456 default: VASSERT(0);
1457 }
1458 return (ir);
1459 }
1460
1461 /*
1462 * ***************************************************************************
1463 * Routine: Aprx_nodeIndexTT
1464 *
1465 * Purpose: Given an index of node DF within the simplex and
1466 * its TT structure, compute global index for that node
1467 *
1468 * Author: Michael Holst, Oleg Korobkin (updated by Ryan Szypowski)
1469 * ***************************************************************************
1470 */
Aprx_nodeIndexTT(Aprx * thee,Re * re,TT * t,int idf)1471 VPUBLIC int Aprx_nodeIndexTT (Aprx *thee, Re *re, TT *t, int idf)
1472 {
1473 int d, thr, ni, nq;
1474
1475 thr = 0;
1476 /* d is dimension of node DF (0=v,1=e,2=f,3=s) */
1477 for (d=0; d<re->dim; d++) {
1478 thr = re->dimIS[d] * re->numDF[d];
1479 if (idf < thr) break;
1480 idf -= thr;
1481 }
1482 ni = idf / re->numDF[d]; /* number of i-simplex */
1483 nq = idf % re->numDF[d]; /* df within this i-simplex */
1484
1485 /* get i-simplex global number */
1486 switch (d) {
1487 case 0: ni = t->vid[ni]; break;
1488 case 1: ni = t->eid[ni]; break;
1489 case 2: ni = t->fid[ni]; break;
1490 case 3: ni = t->sid; break;
1491 default: VASSERT(0);
1492 }
1493
1494 return (Aprx_nodeIndex (thee, re, ni, d, nq));
1495 }
1496
1497 /*
1498 * ***************************************************************************
1499 * Routine: Aprx_deform
1500 *
1501 * Purpose: Deform the mesh.
1502 *
1503 * Author: Michael Holst
1504 * ***************************************************************************
1505 */
Aprx_deform(Aprx * thee,Bvec * def)1506 VPUBLIC int Aprx_deform(Aprx *thee, Bvec *def)
1507 {
1508 int rc,i,numB;
1509 double *defX[MAXV];
1510
1511 /* first make sure "u" makes sense as a mesh displacement */
1512 numB = Aprx_numB(thee);
1513 VWARN( numB == Gem_dimII(thee->gm) );
1514
1515 /* set up the data fields for use by the deformation routine */
1516 for (i=0; i<MAXV; i++) {
1517 defX[i] = VNULL;
1518 }
1519 for (i=0; i<numB; i++) {
1520 defX[i] = Bvec_addrB(def,i);
1521 }
1522
1523 /* call deform */
1524 rc = Gem_deform(thee->gm,numB,defX);
1525
1526 /* return */
1527 return rc;
1528 }
1529
1530