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: am.c,v 1.58 2010/08/12 05:19:12 fetk Exp $"
21 * ***************************************************************************
22 */
23
24 /*
25 * ***************************************************************************
26 * File: am.c
27 *
28 * Purpose: Class AM: methods.
29 *
30 * Author: Michael Holst
31 * ***************************************************************************
32 */
33
34 #include "am_p.h"
35
36 VEMBED(rcsid="$Id: am.c,v 1.58 2010/08/12 05:19:12 fetk Exp $")
37
38 /*
39 * ***************************************************************************
40 * Class AM: Inlineable methods
41 * ***************************************************************************
42 */
43 #if !defined(VINLINE_NAM)
44
45 #endif /* if !defined(VINLINE_NAM) */
46 /*
47 * ***************************************************************************
48 * Class AM: Non-inlineable methods
49 * ***************************************************************************
50 */
51
52 /*
53 * ***************************************************************************
54 * Routine: AM_ctor
55 *
56 * Purpose: The AM constructor.
57 *
58 * Author: Michael Holst
59 * ***************************************************************************
60 */
AM_ctor(Vmem * vmem,Aprx * aprx)61 VPUBLIC AM* AM_ctor(Vmem *vmem, Aprx *aprx)
62 {
63 AM *thee = VNULL;
64
65 VDEBUGIO("AM_ctor: CREATING object..");
66
67 thee = Vmem_malloc( VNULL, 1, sizeof(AM) );
68 if (vmem == VNULL) {
69 thee->vmem = Vmem_ctor( "AM" );
70 thee->iMadeVmem = 1;
71 } else {
72 thee->vmem = vmem;
73 thee->iMadeVmem = 0;
74 }
75
76 /* hand-off of aprx pointer */
77 VASSERT( VNULL != (thee->aprx = aprx) );
78
79 /* no algebraic structures exist yet */
80 thee->algExist = 0;
81
82 /* basic structure data */
83 thee->P = VNULL;
84 thee->A = VNULL;
85 thee->M = VNULL;
86 thee->f = VNULL;
87 thee->u = VNULL;
88 thee->ud = VNULL;
89 thee->ui = VNULL;
90 thee->ut = VNULL;
91 thee->r = VNULL;
92 thee->w0 = VNULL;
93
94 VDEBUGIO("..done.\n");
95
96 return thee;
97 }
98
99 /*
100 * ***************************************************************************
101 * Routine: AM_dtor
102 *
103 * Purpose: The AM destructor.
104 *
105 * Author: Michael Holst
106 * ***************************************************************************
107 */
AM_dtor(AM ** thee)108 VPUBLIC void AM_dtor(AM **thee)
109 {
110 VASSERT( (*thee) != VNULL );
111 if ((*thee) != VNULL) {
112
113 /* destroy the block matrix storage arrays */
114 if ( (*thee)->algExist ) {
115 AM_destroy(*thee);
116 }
117
118 VDEBUGIO("AM_dtor: DESTROYING object..");
119 if ((*thee)->iMadeVmem) Vmem_dtor( &((*thee)->vmem) );
120 Vmem_free( VNULL, 1, sizeof(AM), (void**)thee );
121 VDEBUGIO("..done.\n");
122
123 (*thee) = VNULL;
124 }
125 }
126
127 /*
128 * ***************************************************************************
129 * Routine: AM_create
130 *
131 * Purpose: Create the following internal Alg structures:
132 *
133 * A ==> The tangent matrix (linearization operator)
134 * M ==> The mass matrix
135 * W[] ==> The node vectors
136 *
137 * Author: Michael Holst
138 * ***************************************************************************
139 */
AM_create(AM * thee)140 VPUBLIC void AM_create(AM *thee)
141 {
142 int numB;
143 int numR[MAXV], numO[MAXV][MAXV], *IJA[MAXV][MAXV];
144 MATsym sym[MAXV][MAXV];
145 MATmirror mirror[MAXV][MAXV];
146 MATformat frmt[MAXV][MAXV];
147
148 if ( !thee->algExist ) {
149
150 /* initialization key */
151 thee->algExist = 1;
152
153 /* node interaction; gives structure of matrix */
154 Aprx_interactBlock( thee->aprx, thee->aprx->re, thee->aprx->re,
155 sym, mirror, frmt, numR, numO, IJA );
156
157 /* build sparse (square) block matrices (only integer structures) */
158 numB = Aprx_numB(thee->aprx);
159 thee->A = Bmat_ctor( thee->vmem, "A", numB, numR, numR, mirror );
160 Bmat_initStructure( thee->A, frmt, sym, numO, IJA );
161
162 /* now handle mass matrix (purposely use same mirror struct as A) */
163 thee->M = Bmat_ctor( thee->vmem, "M", numB, numR, numR, mirror );
164 Bmat_copyStructure( thee->M, thee->A );
165
166 /* insert the boundary node information into A and M */
167 Aprx_buildBRC( thee->aprx, thee->A, thee->M );
168
169 /* NODE work vector construction */
170 thee->f = Bvec_ctor( thee->vmem, "f", numB, numR );
171 thee->u = Bvec_ctor( thee->vmem, "u", numB, numR );
172 thee->ud = Bvec_ctor( thee->vmem, "ud", numB, numR );
173 thee->ui = Bvec_ctor( thee->vmem, "ui", numB, numR );
174 thee->ut = Bvec_ctor( thee->vmem, "ut", numB, numR );
175 thee->r = Bvec_ctor( thee->vmem, "r", numB, numR );
176 thee->w0 = Bvec_ctor( thee->vmem, "w0", numB, numR );
177
178 /* get a pointer to prolongation and Bnode objects in geometry object */
179 thee->P = Aprx_P( thee->aprx );
180
181 /* evaluate the trace function u_d (and also u_i and u_t) */
182 Aprx_evalTrace( thee->aprx, thee->ud, thee->ui, thee->ut );
183 }
184 }
185
186 /*
187 * ***************************************************************************
188 * Routine: AM_destroy
189 *
190 * Purpose: Destroy the following internal Alg structures:
191 *
192 * A ==> The tangent matrix (linearization operator)
193 * M ==> The mass matrix
194 * W[] ==> The node vectors
195 *
196 * Author: Michael Holst
197 * ***************************************************************************
198 */
AM_destroy(AM * thee)199 VPUBLIC void AM_destroy(AM *thee)
200 {
201 if (thee->algExist) {
202
203 /* reset the initialization key */
204 thee->algExist = 0;
205
206 /* kill the objects */
207 Bmat_dtor( &(thee->A) );
208 Bmat_dtor( &(thee->M) );
209 Bvec_dtor( &(thee->f) );
210 Bvec_dtor( &(thee->u) );
211 Bvec_dtor( &(thee->ud) );
212 Bvec_dtor( &(thee->ui) );
213 Bvec_dtor( &(thee->ut) );
214 Bvec_dtor( &(thee->r) );
215 Bvec_dtor( &(thee->w0) );
216
217 }
218 }
219
220 /*
221 * ***************************************************************************
222 * Routine: AM_markRefine
223 *
224 * Purpose: Mark a given mesh for refinement.
225 *
226 * Author: Michael Holst
227 * ***************************************************************************
228 */
AM_markRefine(AM * thee,int key,int color,int bkey,double elevel)229 VPUBLIC int AM_markRefine(AM *thee, int key, int color,
230 int bkey, double elevel)
231 {
232 int rc;
233
234 /* make sure we have structure */
235 AM_create(thee);
236
237 /* mark the mesh */
238 rc = Aprx_markRefine( thee->aprx, key, color, bkey, elevel,
239 thee->u, thee->ud, thee->f, thee->r );
240
241 /* return */
242 return rc;
243 }
244
245 /*
246 * ***************************************************************************
247 * Routine: AM_refine
248 *
249 * Purpose: Refine the mesh.
250 *
251 * Author: Michael Holst
252 * ***************************************************************************
253 */
AM_refine(AM * thee,int rkey,int bkey,int pkey)254 VPUBLIC int AM_refine(AM *thee, int rkey, int bkey, int pkey)
255 {
256 int rc;
257
258 /* destroy everything */
259 AM_destroy(thee);
260
261 /* refinement */
262 rc = Aprx_refine(thee->aprx, rkey, bkey, pkey);
263
264 /* recreate everything */
265 AM_create(thee);
266
267 /* return */
268 return rc;
269 }
270
271 /*
272 * ***************************************************************************
273 * Routine: AM_unRefine
274 *
275 * Purpose: Un-refine the mesh.
276 *
277 * Author: Michael Holst
278 * ***************************************************************************
279 */
AM_unRefine(AM * thee,int rkey,int pkey)280 VPUBLIC int AM_unRefine(AM *thee, int rkey, int pkey)
281 {
282 int rc;
283
284 /* destroy everything */
285 AM_destroy(thee);
286
287 /* unrefinement */
288 rc = Aprx_unRefine(thee->aprx, rkey, pkey);
289
290 /* recreate everything */
291 AM_create(thee);
292
293 /* return */
294 return rc;
295 }
296
297 /*
298 * ***************************************************************************
299 * Routine: AM_deform
300 *
301 * Purpose: Deform the mesh.
302 *
303 * Author: Michael Holst
304 * ***************************************************************************
305 */
AM_deform(AM * thee)306 VPUBLIC int AM_deform(AM *thee)
307 {
308 int rc;
309
310 /* deformation */
311 rc = Aprx_deform(thee->aprx, thee->u);
312
313 /* return */
314 return rc;
315 }
316
317 /*
318 * ***************************************************************************
319 * Routine: AM_read
320 *
321 * Purpose: Read in the user-specified initial mesh given in the
322 * "MCSF" or "MCEF" format, and transform into our internal
323 * datastructures.
324 *
325 * Do a little more than a "Aprx_read", in that we also
326 * initialize the extrinsic and intrinsic spatial dimensions
327 * corresponding to the input mesh, and we also then build the
328 * reference elements.
329 *
330 * Notes: See the documentation to Aprx_read for a description of the
331 * mesh input data file format.
332 *
333 * Author: Michael Holst
334 * ***************************************************************************
335 */
AM_read(AM * thee,int key,Vio * sock)336 VPUBLIC int AM_read(AM *thee, int key, Vio *sock)
337 {
338 int rc;
339
340 /* destroy everything */
341 AM_destroy(thee);
342
343 /* read in the mesh; return error if problems */
344 rc = Aprx_read(thee->aprx,key,sock);
345
346 /* recreate everything */
347 AM_create(thee);
348
349 /* return */
350 return rc;
351 }
352
353 /*
354 * ***************************************************************************
355 * Routine: AM_assem
356 *
357 * Purpose: Assemble the linearized problem at a given level.
358 *
359 * Author: Michael Holst
360 * ***************************************************************************
361 */
AM_assem(AM * thee,int evalKey,int energyKey,int residKey,int tangKey,int massKey,int bumpKey,Bvec * u,Bvec * ud,Bvec * f,int ip[],double rp[])362 VPUBLIC double AM_assem(AM *thee,
363 int evalKey, int energyKey, int residKey, int tangKey, int massKey,
364 int bumpKey,
365 Bvec *u, Bvec *ud, Bvec *f, int ip[], double rp[])
366 {
367 double energy = 0.0;
368
369 /* make sure we have structure */
370 AM_create(thee);
371
372 /* assemble the linearized problem */
373 energy = Aprx_assem( thee->aprx,
374 evalKey, energyKey, residKey, tangKey, massKey,
375 bumpKey,
376 ip, rp, thee->A, thee->M, f, u, ud );
377
378 /* return */
379 return energy;
380 }
381
382 /*
383 * ***************************************************************************
384 * Routine: AM_evalJ
385 *
386 * Purpose: Assemble the energy functional at the current solution.
387 *
388 * Author: Michael Holst
389 * ***************************************************************************
390 */
AM_evalJ(AM * thee)391 VPUBLIC double AM_evalJ(AM *thee)
392 {
393 int ip[10], evalKey, energyKey, residKey, tangKey, massKey, bumpKey;
394 double rp[10], energy;
395
396 /* make sure we have structure */
397 AM_create(thee);
398
399 /* compute the energy */
400 evalKey = 1;
401 energyKey = 1;
402 residKey = 0;
403 tangKey = 0;
404 massKey = 0;
405 bumpKey = 0;
406 energy = AM_assem(thee,
407 evalKey, energyKey, residKey, tangKey, massKey,
408 bumpKey,
409 thee->u, thee->ud, thee->f, ip, rp);
410
411 /* return */
412 return energy;
413 }
414
415 /*
416 * ***************************************************************************
417 * Routine: AM_evalFunc
418 *
419 * Purpose: Evaluate a finite element function.
420 *
421 * Author: Michael Holst
422 * ***************************************************************************
423 */
AM_evalFunc(AM * thee,int number,int block,int numPts,double * pts,double * vals,int * marks)424 VPUBLIC void AM_evalFunc(AM *thee,
425 int number, int block, int numPts, double *pts,
426 double *vals, int *marks)
427 {
428 /* make sure we have structure */
429 AM_create(thee);
430
431 /* function evaluation */
432 Aprx_evalFunc( thee->aprx, thee->u, block, numPts, pts, vals, marks );
433 }
434
435 /*
436 * ***************************************************************************
437 * Routine: AM_bndIntegral
438 *
439 * Purpose: Perform a boundary integral.
440 *
441 * Author: Michael Holst
442 * ***************************************************************************
443 */
AM_bndIntegral(AM * thee)444 VPUBLIC void AM_bndIntegral(AM *thee)
445 {
446 /* make sure we have structure */
447 AM_create(thee);
448
449 /* boundary integral evaluation */
450 Aprx_bndIntegral( thee->aprx, thee->u, thee->ud, thee->ut );
451 }
452
453 /*
454 * ***************************************************************************
455 * Routine: AM_evalError
456 *
457 * Purpose: Evaluate error in the current solution.
458 *
459 * Author: Michael Holst
460 * ***************************************************************************
461 */
AM_evalError(AM * thee,int pcolor,int key)462 VPUBLIC double AM_evalError(AM *thee, int pcolor, int key)
463 {
464 int rc;
465
466 /* make sure we have structure */
467 AM_create(thee);
468
469 /* error evaluation */
470 rc = Aprx_evalError( thee->aprx,
471 pcolor, key, thee->u, thee->ud, thee->ut );
472
473 /* return */
474 return rc;
475 }
476
477 /*
478 * ***************************************************************************
479 * Routine: AM_applyDiriZero
480 *
481 * Purpose: Apply zero dirichlet condition at a given level.
482 *
483 * Author: Michael Holst
484 * ***************************************************************************
485 */
AM_applyDiriZero(AM * thee,Bvec * v)486 VPUBLIC void AM_applyDiriZero(AM *thee, Bvec *v)
487 {
488 /* make sure we have structure */
489 AM_create(thee);
490
491 /* apply boundary condition */
492 Bvec_bnd( v, thee->A, 0 );
493 }
494
495 /*
496 * ***************************************************************************
497 * Routine: AM_iniGuess
498 *
499 * Purpose: Setup an initial guess at a given level.
500 *
501 * Author: Michael Holst
502 * ***************************************************************************
503 */
AM_iniGuess(AM * thee,Bvec * v)504 VPUBLIC void AM_iniGuess(AM *thee, Bvec *v)
505 {
506 int i,j;
507
508 /* make sure we have structure */
509 AM_create(thee);
510
511 /* apply initial condition */
512 for (i=0; i<Bvec_numB(v); i++)
513 for (j=0; j<Bvec_numRB(v,i); j++)
514 Bvec_setB( v,i,j,Bvec_valB(thee->ui,i,j) );
515 }
516
517 /*
518 * ***************************************************************************
519 * Routine: AM_part
520 *
521 * Purpose: Partition the mesh using the matching Alg level.
522 *
523 * Author: Michael Holst
524 * ***************************************************************************
525 */
AM_part(AM * thee,int pkey,int pwht,int ppow)526 VPUBLIC int AM_part(AM *thee, int pkey, int pwht, int ppow)
527 {
528 /* make sure we have structure */
529 AM_create(thee);
530
531 /* partition the mesh */
532 return Aprx_part( thee->aprx, pkey, pwht, ppow );
533 }
534
535 /*
536 * ***************************************************************************
537 * Routine: AM_partSet
538 *
539 * Purpose: Set the partition color.
540 *
541 * Author: Michael Holst
542 * ***************************************************************************
543 */
AM_partSet(AM * thee,int pcolor)544 VPUBLIC int AM_partSet(AM *thee, int pcolor)
545 {
546 /* make sure we have structure */
547 AM_create(thee);
548
549 /* set partition */
550 return Aprx_partSet( thee->aprx, pcolor );
551 }
552
553 /*
554 * ***************************************************************************
555 * Routine: AM_partSmooth
556 *
557 * Purpose: Do a partition smoothing.
558 *
559 * Author: Michael Holst
560 * ***************************************************************************
561 */
AM_partSmooth(AM * thee)562 VPUBLIC int AM_partSmooth(AM *thee)
563 {
564 /* make sure we have structure */
565 AM_create(thee);
566
567 /* smooth partition */
568 return Aprx_partSmooth( thee->aprx );
569 }
570
571 /*
572 * ***************************************************************************
573 * Routine: AM_printJ
574 *
575 * Purpose: Print the energy.
576 *
577 * Author: Michael Holst
578 * ***************************************************************************
579 */
AM_printJ(AM * thee)580 VPUBLIC void AM_printJ(AM *thee)
581 {
582 double energy;
583
584 /* make sure we have structure */
585 AM_create(thee);
586
587 /* print energy */
588 energy = AM_evalJ(thee);
589 Vnm_print(2,"AM_printJ: Energy=<%g>\n", energy);
590 }
591
592 /*
593 * ***************************************************************************
594 * Routine: AM_printA
595 *
596 * Purpose: Print the system matrix.
597 *
598 * Author: Michael Holst
599 * ***************************************************************************
600 */
AM_printA(AM * thee)601 VPUBLIC void AM_printA(AM *thee)
602 {
603 /* make sure we have structure */
604 AM_create(thee);
605
606 /* print A */
607 Bmat_print( thee->A );
608 }
609
610 /*
611 * ***************************************************************************
612 * Routine: AM_printAnoD
613 *
614 * Purpose: Print the system matrix with Dirichlet rows/cols zeroed.
615 *
616 * Author: Michael Holst
617 * ***************************************************************************
618 */
AM_printAnoD(AM * thee)619 VPUBLIC void AM_printAnoD(AM *thee)
620 {
621 /* make sure we have structure */
622 AM_create(thee);
623
624 /* print A */
625 Bmat_printNoD( thee->A );
626 }
627
628 /*
629 * ***************************************************************************
630 * Routine: AM_printAsp
631 *
632 * Purpose: Print the system matrix in MATLAB sparse format.
633 *
634 * Author: Michael Holst
635 * ***************************************************************************
636 */
AM_printAsp(AM * thee,char * fname)637 VPUBLIC void AM_printAsp(AM *thee, char *fname)
638 {
639 /* make sure we have structure */
640 AM_create(thee);
641
642 /* print A */
643 Bmat_printSp( thee->A, fname, 0 );
644 }
645
646 /*
647 * ***************************************************************************
648 * Routine: AM_printAspNoD
649 *
650 * Purpose: Print the system matrix in MATLAB sparse format with
651 * Dirichlet rows/cols zeroed.
652 *
653 * Author: Michael Holst
654 * ***************************************************************************
655 */
AM_printAspNoD(AM * thee,char * fname)656 VPUBLIC void AM_printAspNoD(AM *thee, char *fname)
657 {
658 /* make sure we have structure */
659 AM_create(thee);
660
661 /* print A */
662 Bmat_printSpNoD( thee->A, fname, 0 );
663 }
664
665 /*
666 * ***************************************************************************
667 * Routine: AM_printP
668 *
669 * Purpose: Print the prolongation matrix.
670 *
671 * Author: Michael Holst
672 * ***************************************************************************
673 */
AM_printP(AM * thee)674 VPUBLIC void AM_printP(AM *thee)
675 {
676 /* make sure we have structure */
677 AM_create(thee);
678
679 /* print P */
680 Bmat_print( thee->P );
681 }
682
683 /*
684 * ***************************************************************************
685 * Routine: AM_printPsp
686 *
687 * Purpose: Print the prolongation matrix in MATLAB sparse format.
688 *
689 * Author: Michael Holst
690 * ***************************************************************************
691 */
AM_printPsp(AM * thee,char * fname)692 VPUBLIC void AM_printPsp(AM *thee, char *fname)
693 {
694 /* make sure we have structure */
695 AM_create(thee);
696
697 /* print P */
698 Bmat_printSp( thee->P, fname, 0 );
699 }
700
701 /*
702 * ***************************************************************************
703 * Routine: AM_printV
704 *
705 * Purpose: Print a specified vector.
706 *
707 * Author: Michael Holst
708 * ***************************************************************************
709 */
AM_printV(AM * thee,int number)710 VPUBLIC void AM_printV(AM *thee, int number)
711 {
712 int maxw, W_f, W_u, W_ud, W_ui, W_ut, W_r, W_w0;
713 Bvec *W[7];
714 char *W_name[7];
715
716 /* make sure we have structure */
717 AM_create(thee);
718
719 maxw = 7;
720 W_f = 0, /* rhs/residual */
721 W_u = 1, /* solution/error */
722 W_ud = 2, /* dirichlet function */
723 W_ui = 3, /* interior dirichlet function (initial guess) */
724 W_ut = 4, /* analytical solution (for testing) */
725 W_r = 5, /* residual/temp vector */
726 W_w0 = 6, /* work vector */
727
728 W[W_f] = thee->f;
729 W[W_u] = thee->u;
730 W[W_ud] = thee->ud;
731 W[W_ui] = thee->ui;
732 W[W_ut] = thee->ut;
733 W[W_r] = thee->r;
734 W[W_w0] = thee->w0;
735
736 W_name[W_f] = "f";
737 W_name[W_u] = "u";
738 W_name[W_ud] = "ud";
739 W_name[W_ui] = "ui";
740 W_name[W_ut] = "ut";
741 W_name[W_r] = "r";
742 W_name[W_w0] = "w0";
743
744 if (number == -W_ud) number = W_ud;
745 if ((-maxw < number) && (number < 0)) {
746 Vnm_print(0,"AM_printV: forming <w0=%s+ud> to print <w0>\n",
747 W_name[-number]);
748 Bvec_copy( W[W_w0], W[-number] );
749 Bvec_axpy( W[W_w0], W[W_ud], 1. );
750 Bvec_print( W[W_w0] );
751 } else if ((0 <= number) && (number < maxw)) {
752 Vnm_print(0,"AM_printV: printing <%s>\n", W_name[number]);
753 Bvec_print( W[number] );
754 }
755 }
756
757 /*
758 * ***************************************************************************
759 * Routine: AM_printVsp
760 *
761 * Purpose: Print a vector in MATLAB sparse format.
762 *
763 * Author: Michael Holst
764 * ***************************************************************************
765 */
AM_printVsp(AM * thee,int number,char * fname)766 VPUBLIC void AM_printVsp(AM *thee, int number, char *fname)
767 {
768 int maxw, W_f, W_u, W_ud, W_ui, W_ut, W_r, W_w0;
769 Bvec *W[7];
770 char *W_name[7];
771
772 /* make sure we have structure */
773 AM_create(thee);
774
775 maxw = 7;
776 W_f = 0, /* rhs/residual */
777 W_u = 1, /* solution/error */
778 W_ud = 2, /* dirichlet function */
779 W_ui = 3, /* interior dirichlet function (initial guess) */
780 W_ut = 4, /* analytical solution (for testing) */
781 W_r = 5, /* residual/temp vector */
782 W_w0 = 6, /* work vector */
783
784 W[W_f] = thee->f;
785 W[W_u] = thee->u;
786 W[W_ud] = thee->ud;
787 W[W_ui] = thee->ui;
788 W[W_ut] = thee->ut;
789 W[W_r] = thee->r;
790 W[W_w0] = thee->w0;
791
792 W_name[W_f] = "f";
793 W_name[W_u] = "u";
794 W_name[W_ud] = "ud";
795 W_name[W_ui] = "ui";
796 W_name[W_ut] = "ut";
797 W_name[W_r] = "r";
798 W_name[W_w0] = "w0";
799
800 if (number == -W_ud) number = W_ud;
801 if ((-maxw < number) && (number < 0)) {
802 Vnm_print(0,"AM_printVsp: forming <w0=%s+ud> to print <w0>\n",
803 W_name[-number]);
804 Bvec_copy( W[W_w0], W[-number] );
805 Bvec_axpy( W[W_w0], W[W_ud], 1. );
806 Bvec_printSp( W[W_w0], fname );
807 } else if ((0 <= number) && (number < maxw)) {
808 Vnm_print(0,"AM_printVsp: printing <%s>\n", W_name[number]);
809 Bvec_printSp( W[number], fname );
810 }
811 }
812
813 /*
814 * ***************************************************************************
815 * Routine: AM_writeGEOM
816 *
817 * Purpose: Write out a mesh in some format.
818 *
819 * Author: Michael Holst
820 * ***************************************************************************
821 */
AM_writeGEOM(AM * thee,Vio * sock,int defKey,int colKey,int chartKey,double gluVal,int fkey,int number,char * format)822 VPUBLIC void AM_writeGEOM(AM *thee, Vio *sock,
823 int defKey, int colKey, int chartKey, double gluVal, int fkey,
824 int number, char *format)
825 {
826 int maxw, W_f, W_u, W_ud, W_ui, W_ut, W_r, W_w0;
827 Bvec *W[7];
828 char *W_name[7];
829
830 /* make sure we have structure */
831 AM_create(thee);
832
833 maxw = 7;
834 W_f = 0, /* rhs/residual */
835 W_u = 1, /* solution/error */
836 W_ud = 2, /* dirichlet function */
837 W_ui = 3, /* interior dirichlet function (initial guess) */
838 W_ut = 4, /* analytical solution (for testing) */
839 W_r = 5, /* residual/temp vector */
840 W_w0 = 6, /* work vector */
841
842 W[W_f] = thee->f;
843 W[W_u] = thee->u;
844 W[W_ud] = thee->ud;
845 W[W_ui] = thee->ui;
846 W[W_ut] = thee->ut;
847 W[W_r] = thee->r;
848 W[W_w0] = thee->w0;
849
850 W_name[W_f] = "f";
851 W_name[W_u] = "u";
852 W_name[W_ud] = "ud";
853 W_name[W_ui] = "ui";
854 W_name[W_ut] = "ut";
855 W_name[W_r] = "r";
856 W_name[W_w0] = "w0";
857
858 if (number == -W_ud) number = W_ud;
859 if ((-maxw < number) && (number < 0)) {
860 Vnm_print(0,"AM_writeGEOM: forming <w0=%s+ud> to print <w0>\n",
861 W_name[-number]);
862 Bvec_copy( W[W_w0], W[-number] );
863 Bvec_axpy( W[W_w0], W[W_ud], 1. );
864 Aprx_writeGEOM( thee->aprx,
865 sock, defKey, colKey, chartKey, gluVal, fkey, W[W_w0], format );
866 } else if ((0 <= number) && (number < maxw)) {
867 Vnm_print(0,"AM_writeGEOM: printing <%s>\n", W_name[number]);
868 Aprx_writeGEOM( thee->aprx,
869 sock, defKey, colKey, chartKey, gluVal, fkey, W[number], format );
870 }
871 }
872
873 /*
874 * ***************************************************************************
875 * Routine: AM_writeSOL
876 *
877 * Purpose: Write out a solution in some format.
878 *
879 * Author: Michael Holst
880 * ***************************************************************************
881 */
AM_writeSOL(AM * thee,Vio * sock,int number,char * format)882 VPUBLIC void AM_writeSOL(AM *thee, Vio *sock, int number, char *format)
883 {
884 int maxw, W_f, W_u, W_ud, W_ui, W_ut, W_r, W_w0;
885 Bvec *W[7];
886 char *W_name[7];
887
888 /* make sure we have structure */
889 AM_create(thee);
890
891 maxw = 7;
892 W_f = 0, /* rhs/residual */
893 W_u = 1, /* solution/error */
894 W_ud = 2, /* dirichlet function */
895 W_ui = 3, /* interior dirichlet function (initial guess) */
896 W_ut = 4, /* analytical solution (for testing) */
897 W_r = 5, /* residual/temp vector */
898 W_w0 = 6, /* work vector */
899
900 W[W_f] = thee->f;
901 W[W_u] = thee->u;
902 W[W_ud] = thee->ud;
903 W[W_ui] = thee->ui;
904 W[W_ut] = thee->ut;
905 W[W_r] = thee->r;
906 W[W_w0] = thee->w0;
907
908 W_name[W_f] = "f";
909 W_name[W_u] = "u";
910 W_name[W_ud] = "ud";
911 W_name[W_ui] = "ui";
912 W_name[W_ut] = "ut";
913 W_name[W_r] = "r";
914 W_name[W_w0] = "w0";
915
916 if (number == -W_ud) number = W_ud;
917 if ((-maxw < number) && (number < 0)) {
918 Vnm_print(0,"AM_writeSOL: forming <w0=%s+ud> to print <w0>\n",
919 W_name[-number]);
920 Bvec_copy( W[W_w0], W[-number] );
921 Bvec_axpy( W[W_w0], W[W_ud], 1. );
922 Aprx_writeSOL( thee->aprx, sock, W[W_w0], format );
923 } else if ((0 <= number) && (number < maxw)) {
924 Vnm_print(0,"AM_writeSOL: printing <%s>\n", W_name[number]);
925 Aprx_writeSOL( thee->aprx, sock, W[number], format );
926 }
927 }
928
929 /*
930 * ***************************************************************************
931 * Routine: AM_memChk
932 *
933 * Purpose: Print the exact current malloc usage.
934 *
935 * Author: Michael Holst
936 * ***************************************************************************
937 */
AM_memChk(AM * thee)938 VPUBLIC void AM_memChk(AM *thee)
939 {
940 if (thee->iMadeVmem) Vmem_print(thee->vmem);
941 if (thee->P != VNULL) Bmat_memChk(thee->P);
942 if (thee->A != VNULL) Bmat_memChk(thee->A);
943 if (thee->M != VNULL) Bmat_memChk(thee->M);
944 if (thee->f != VNULL) Bvec_memChk(thee->f);
945 if (thee->u != VNULL) Bvec_memChk(thee->u);
946 if (thee->ud != VNULL) Bvec_memChk(thee->ud);
947 if (thee->ui != VNULL) Bvec_memChk(thee->ui);
948 if (thee->ut != VNULL) Bvec_memChk(thee->ut);
949 if (thee->r != VNULL) Bvec_memChk(thee->r);
950 if (thee->w0 != VNULL) Bvec_memChk(thee->w0);
951 }
952
953