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: assem.c,v 1.21 2010/08/12 05:18:16 fetk Exp $"
21 * ***************************************************************************
22 */
23
24 /*
25 * ***************************************************************************
26 * File: assem.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: assem.c,v 1.21 2010/08/12 05:18:16 fetk Exp $")
37
38 /*
39 * ***************************************************************************
40 * Class Alg: Inlineable methods
41 * ***************************************************************************
42 */
43 #if !defined(VINLINE_APRX)
44
45 #endif /* if !defined(VINLINE_APRX) */
46 /*
47 * ***************************************************************************
48 * Class Alg: Non-inlineable methods
49 * ***************************************************************************
50 */
51
52 /*
53 * ***************************************************************************
54 * Routine: Aprx_buildXq
55 *
56 * Purpose: Translate reference element coordinate to the current triangle.
57 *
58 * Author: Michael Holst and Ryan Szypowski
59 * ***************************************************************************
60 */
Aprx_buildXq(Aprx * thee,Re * re,int qp,int face,TT * t,double xq[])61 VPUBLIC void Aprx_buildXq(Aprx *thee, Re *re,
62 int qp, int face, TT *t,
63 double xq[])
64 {
65 int i, j;
66
67 Gem *gm;
68 gm = thee->gm; /* shortcut to gem structure */
69
70 /* Get quad pt by mapping master el quad pt to this el */
71 for (i=0; i<Gem_dimII(gm); i++) {
72 xq[i] = t->bb[i];
73 for (j=0;j<Gem_dimII(gm);j++)
74 xq[i] += ( t->ff[i][j] * Re_x_hi(re,qp,j,face) );
75 }
76 }
77
78 /*
79 * ***************************************************************************
80 * Routine: Aprx_buildPhi
81 *
82 * Purpose: Build the scalar basis functions and their gradients
83 * on this element.
84 *
85 * Author: Michael Holst and Ryan Szypowski
86 * ***************************************************************************
87 */
Aprx_buildPhi(Aprx * thee,Re * re,int qp,int face,TT * t,double phi[],double phix[][3])88 VPUBLIC void Aprx_buildPhi(Aprx *thee, Re *re,
89 int qp, int face, TT *t,
90 double phi[], double phix[][3])
91 {
92 int i, j, k;
93
94 Gem *gm;
95 gm = thee->gm; /* shortcut to gem structure */
96
97 /* Get basis functions; transform grads to arbitrary elm */
98 for (i=0; i<Re_numP(re,-1); i++) {
99 phi[i] = Re_phi_hi(re,qp,i,face);
100 t->phi[i] = phi[i];
101 for (j=0; j<Gem_dim(gm); j++) {
102 phix[i][j] = 0.;
103 for (k=0; k<Gem_dim(gm); k++)
104 phix[i][j] += ( t->gg[k][j] * Re_phix2_hi(re,qp,i,k,face) );
105 t->phix[i][j] = phix[i][j];
106 }
107 }
108 }
109
110 /*
111 * ***************************************************************************
112 * Routine: Aprx_buildU
113 *
114 * Purpose: Build an interpolant of the solution and its gradient.
115 *
116 * Notes: evalKey==0 ==> use U=[0+ud] as evaluation point
117 * evalKey>=1 ==> use U=[u+ud] as evaluation point
118 *
119 * Author: Michael Holst and Ryan Szypowski
120 * ***************************************************************************
121 */
Aprx_buildU(Aprx * thee,Re * re,int evalKey,int qp,int face,TT * t,double U[],double dU[][3],Bvec * Wu,Bvec * Wud)122 VPUBLIC void Aprx_buildU(Aprx *thee, Re *re,
123 int evalKey, int qp, int face, TT *t,
124 double U[], double dU[][3],
125 Bvec *Wu, Bvec *Wud)
126 {
127 int i, j, k, gi;
128 double u_u[VMAXP][MAXV], u_ud[VMAXP][MAXV];
129 double phi[VMAXP], phix[VMAXP][3];
130
131 Gem *gm;
132 gm = thee->gm; /* shortcut to gem structure */
133
134 /* Get basis functions and gradients */
135 Aprx_buildPhi(thee, re, qp, face, t, phi, phix);
136
137 /* Handle linearization differences between zero and non-zero functions */
138 if (evalKey == 0) {
139 for (j=0; j<Re_numP(re,-1); j++) {
140 gi = Aprx_nodeIndexTT (thee, re, t, j);
141 for (i=0; i<PDE_vec(thee->pde); i++) {
142 u_u[j][i] = 0.;
143 u_ud[j][i] = Bvec_valB( Wud, i, gi );
144 }
145 }
146 } else {
147 for (j=0; j<Re_numP(re,-1); j++) {
148 gi = Aprx_nodeIndexTT (thee, re, t, j);
149 for (i=0; i<PDE_vec(thee->pde); i++) {
150 u_u[j][i] = Bvec_valB( Wu, i, gi );
151 u_ud[j][i] = Bvec_valB( Wud, i, gi );
152 }
153 }
154 }
155
156 /* Initialize function [U+UD] and its gradient [dU+dUD] */
157 for (i=0; i<PDE_vec(thee->pde); i++) {
158 U[i] = 0.;
159 for (k=0; k<Gem_dim(gm); k++) {
160 dU[i][k] = 0.;
161 for (j=0; j<Re_numP(re,-1); j++) {
162 if (k==0) U[i] += phi[j] * ( u_u[j][i] + u_ud[j][i] );
163 dU[i][k] += phix[j][k] * ( u_u[j][i] + u_ud[j][i] );
164 }
165 }
166 }
167 }
168
169 /*
170 * ***************************************************************************
171 * Routine: Aprx_initSpace
172 *
173 * Purpose: Initialize vector test function "i" and its gradient,
174 * generated by scalar test function "r" for this element.
175 *
176 * Author: Michael Holst
177 * ***************************************************************************
178 */
Aprx_initSpace(Aprx * thee,int i,int r,double phi[],double phix[][3],double V[],double dV[][3])179 VPUBLIC void Aprx_initSpace(Aprx *thee, int i, int r,
180 double phi[], double phix[][3], double V[], double dV[][3])
181 {
182 int j, k;
183 for (j=0; j<PDE_vec(thee->pde); j++) {
184 V[j] = ( (j==i) ? phi[r] : 0. );
185 for (k=0; k<Gem_dim(thee->gm); k++)
186 dV[j][k] = ( (j==i) ? phix[r][k] : 0. );
187 }
188 }
189
190 /*
191 * ***************************************************************************
192 * Routine: Aprx_initEmat
193 *
194 * Purpose: Initialize an element matrix.
195 *
196 * Author: Michael Holst
197 * ***************************************************************************
198 */
Aprx_initEmat(Aprx * thee,int bumpKey,TT * t,Emat * em)199 VPUBLIC void Aprx_initEmat(Aprx *thee, int bumpKey, TT *t, Emat *em)
200 {
201 int r, s, iv, iw, gr;
202 Re *re;
203 Gem *gm;
204
205 if(bumpKey) re = thee->reB[0];
206 else re = thee->re[0];
207 gm = thee->gm;
208
209 em->J = 0.;
210 for (iv=0; iv<PDE_vec(thee->pde); iv++) {
211 for (r=0; r<Re_numP(re,-1); r++) {
212 gr = Aprx_nodeIndexTT (thee, re, t, r);
213 em->NI[iv][r] = gr;
214 #if 0
215 fprintf(stderr,"*** gr=<%d> ***\n", gr);
216 #else
217 if (bumpKey) em->TP[iv][r] = SS_faceType(t->s, r);
218 else em->TP[iv][r] = Bnode_data(thee->B,iv)[gr].type;
219 #endif
220 em->F[iv][r] = 0.;
221 for (iw=0; iw<PDE_vec(thee->pde); iw++) {
222 for (s=0; s<Re_numP(re,-1); s++) {
223 em->A[iv][iw][r][s] = 0.;
224 em->M[iv][iw][r][s] = 0.;
225 }
226 }
227 }
228 }
229 }
230
231 /*
232 * ***************************************************************************
233 * Routine: Aprx_quadEmat
234 *
235 * Purpose: Build an element matrix and an element load vector for a
236 * single given element, using quadrature.
237 *
238 * Notes: We handle both the volume and surface cases here, using
239 * "face" and "t" as follows:
240 *
241 * face < 0 ==> do volume quadrature; we take the node numbers
242 * as just the identity mapping.
243 *
244 * face >= 0 ==> do face quadrature for face "face", where the node
245 * numbers for face are taken from t->loc[face][].
246 *
247 * Author: Michael Holst
248 * ***************************************************************************
249 */
Aprx_quadEmat(Aprx * thee,int evalKey,int energyKey,int residKey,int tangKey,int massKey,int bumpKey,int face,TT * t,Emat * em,Bvec * Wu,Bvec * Wud)250 VPUBLIC void Aprx_quadEmat(Aprx *thee,
251 int evalKey, int energyKey, int residKey, int tangKey, int massKey,
252 int bumpKey,
253 int face, TT *t, Emat *em,
254 Bvec *Wu, Bvec *Wud)
255 {
256 /* some variables */
257 int i, qp, r, s, r2, s2, gr, gs, iv, iw, rrange, sym, key, iw_formed;
258 int loc[VMAXP];
259 double Dw, xq[3], phi[VMAXP], phix[VMAXP][3], gamma, theta, energy, jacD;
260 double U[MAXV], dU[MAXV][3];
261 double W[MAXV], dW[MAXV][3];
262 double V[MAXV], dV[MAXV][3];
263 Re *re, *reU;
264 Gem *gm;
265
266 /* hard-coded element for now... */
267 if(bumpKey) {
268 re = thee->reB[0];
269 if(bumpKey == 2) reU = thee->reB[0];
270 else reU = thee->re[0];
271 }
272 else {
273 re = thee->re[0];
274 reU = thee->re[0];
275 }
276 gm = thee->gm;
277
278 /* deal with differences between a volume and a surface quadrature */
279 if (face<0) {
280 if ((evalKey == 0) || (evalKey == 1)) {
281 key = 0; /* build primal problem: volume form */
282 } else { /* if ((evalKey == 2) || (evalKey == 3)) */
283 key = 2; /* build dual problem: volume form */
284 }
285 jacD = t->D;
286 } else {
287 if ((evalKey == 0) || (evalKey == 1)) {
288 key = 1; /* build primal problem: surface form */
289 } else { /* if ((evalKey == 2) || (evalKey == 3)) */
290 key = 3; /* build dual problem: surface form */
291 }
292 jacD = t->faceD[face];
293 }
294 rrange = Re_numP(re,-1);
295 for (i=0;i<rrange;i++) loc[i] = i;
296
297 /* Go thru quad pts m=1:Qx, computing contrib. to ele matrix. */
298 for (qp=0; qp<Re_numQ_hi(re,face); qp++) {
299
300 /* jacobian and quadrature weight */
301 Dw = jacD * Re_w_hi(re,qp,face);
302
303 /* build basis functions and translate quadrature point to global domain */
304 Aprx_buildPhi(thee,re,qp,face,t,phi,phix);
305 Aprx_buildU(thee,reU,evalKey,qp,face,t,U,dU,Wu,Wud);
306 Aprx_buildXq(thee,re,qp,face,t,xq);
307
308 /* once-per-point PDE initialization */
309 thee->pde->initPoint(thee->pde,key,t->gchart,xq,U,dU);
310
311 /* assemble energy if required */
312 if (energyKey > 0) {
313
314 /* eval (interior or boundary) energy functional */
315 energy = thee->pde->Ju(thee->pde,key);
316
317 /* add in contribution to vector */
318 em->J += Dw*energy;
319 }
320
321 /* outer basis loop; only need if building a residual/load vector */
322 if ((residKey > 0) || (tangKey > 0) || (massKey > 0)) {
323
324 /* begin dealing with inner-products of the form (r,s) */
325 for (r2=0; r2<rrange; r2++) {
326 r = loc[r2];
327 gr = Aprx_nodeIndexTT (thee, re, t, r);
328
329 /* which V-basis function */
330 for (iv=0; iv<PDE_vec(thee->pde); iv++) {
331 t->focusIV = iv;
332
333 /* initialize function/gradient */
334 Aprx_initSpace(thee,iv,r,phi,phix,V,dV);
335
336 /* assemble residual if required */
337 if (residKey > 0) {
338
339 /* eval (interior or boundary) nonlinear form */
340 gamma = thee->pde->Fu_v(thee->pde,key,V,dV);
341
342 /* add in contribution to vector */
343 em->F[iv][r] -= Dw*gamma;
344 }
345
346 /* inner basis loop; only need if building a matrix */
347 if ((tangKey > 0) || (massKey > 0)) {
348
349 /* build only upper-triangle if possible */
350 for (s2=0; s2<rrange; s2++) {
351 s = loc[s2];
352 gs = Aprx_nodeIndexTT (thee, re, t, s);
353
354 /* which W-basis function */
355 for (iw=0; iw<PDE_vec(thee->pde); iw++) {
356 t->focusIW = iw;
357 iw_formed = 0;
358 sym = thee->pde->sym[iv][iw];
359
360 /* assemble tangent matrix if required */
361 if (tangKey > 0) {
362
363 /* handle symmetries in tangent form */
364 if ((sym==0)||((sym==1) && (gr<=gs))) {
365
366 /* initialize function/gradient */
367 Aprx_initSpace(thee,iw,s,phi,phix,W,dW);
368 iw_formed = 1;
369
370 /* eval (interior or bndry) form */
371 theta = thee->pde->DFu_wv(thee->pde,
372 key,W,dW,V,dV);
373
374 /* add in contribution to matrix */
375 em->A[iv][iw][r][s] += Dw*theta;
376
377 } /* symmetries in the tengent form */
378 } /* assemble tangent matrix */
379
380 /* assemble mass matrix if required */
381 if (massKey > 0) {
382
383 /* handle symmetries in the mass form */
384 if ((sym==0)||((sym==1) && (gr<=gs))) {
385
386 /* initialize function/gradient */
387 if (!iw_formed) {
388 Aprx_initSpace(thee,
389 iw,s,phi,phix,W,dW);
390 iw_formed = 1;
391 }
392
393 /* eval (interior or bndry) form */
394 theta = thee->pde->p_wv(thee->pde,
395 key,W,V);
396
397 /* add in contribution to matrix */
398 em->M[iv][iw][r][s] += Dw*theta;
399
400 } /* symmetries in the mass form */
401 } /* assemble mass matrix */
402
403 } /* which W-basis function */
404 } /* s; loop over element matrix columns */
405 } /* if tangent or mass matrix needed */
406 } /* which V-basis function */
407 } /* r; loop over element matrix rows (and, faces of simplex) */
408 } /* if residual/load vector needed */
409 } /* m; loop over quadrature points */
410 }
411
412 /*
413 * ***************************************************************************
414 * Routine: Aprx_fanEmat
415 *
416 * Purpose: Fan out an element matrix to the global matrix,
417 * and fan out an element vector to the global vector.
418 *
419 * Author: Michael Holst
420 * ***************************************************************************
421 */
Aprx_fanEmat(Aprx * thee,int evalKey,int energyKey,int residKey,int tangKey,int massKey,int bumpKey,Emat * em,Bmat * A,Bmat * M,Bvec * F)422 VPUBLIC void Aprx_fanEmat(Aprx *thee,
423 int evalKey, int energyKey, int residKey, int tangKey, int massKey,
424 int bumpKey,
425 Emat *em,
426 Bmat *A, Bmat *M, Bvec *F)
427 {
428 int iv, iw, sym, r, s;
429 #if 0
430 int j;
431 #endif
432 Re *re = thee->re[0];
433
434 if(bumpKey) re = thee->reB[0];
435
436 if ((residKey > 0) || (tangKey > 0) || (massKey > 0)) {
437 for (r=0; r<Re_numP(re,-1); r++) { /* for (r=0; r<Gem_dimVV(thee->gm); r++) */
438 for (iv=0; iv<PDE_vec(thee->pde); iv++) {
439
440 if (residKey > 0) {
441 if (VDIRICHLET( em->TP[iv][r] )) {
442 /* vr is a dirichlet point */
443 Bvec_setB(F, iv, em->NI[iv][r], 0.);
444 } else {
445 /* vr is not a dirichlet point */
446 Bvec_addToB(F, iv, em->NI[iv][r], em->F[iv][r]);
447 }
448 }
449
450 if ((tangKey > 0) || (massKey > 0)) {
451 for (s=0; s<Re_numP(re,-1); s++) { /* for (s=0; s<Gem_dimVV(thee->gm); s++) */
452 for (iw=0; iw<PDE_vec(thee->pde); iw++) {
453 sym = thee->pde->sym[iv][iw];
454
455 /* assemble tangent matrix */
456 if (tangKey > 0) {
457 /* handle symmetries in tangent form */
458 if ( (sym==0) || ((sym==1)
459 &&(em->NI[iv][r]<=em->NI[iw][s])) ) {
460
461 /* one or both vr|vs are dirichlet */
462 if ( VDIRICHLET( em->TP[iv][r] )
463 || VDIRICHLET( em->TP[iw][s] ) ) {
464
465 /* vr is dirichlet */
466 if (VDIRICHLET( em->TP[iv][r] ) ) {
467 #if 1
468 Bmat_set( A,iv,iv, em->NI[iv][r], em->NI[iv][r], 1. );
469 #else
470 for (j=0; j<PDE_vec(thee->pde); j++)
471 Bmat_set( A,j,j,
472 em->NI[iv][r],
473 em->NI[iv][r], 1. );
474 #endif
475 }
476 /* vs is dirichlet */
477 if (VDIRICHLET( em->TP[iw][s] ) ) {
478 #if 1
479 Bmat_set( A,iw,iw, em->NI[iw][s], em->NI[iw][s], 1. );
480 #else
481 for (j=0; j<PDE_vec(thee->pde); j++)
482 Bmat_set( A,j,j,
483 em->NI[iw][s],
484 em->NI[iw][s], 1. );
485 #endif
486 }
487
488 /* neither vr|vs are dirichlet */
489 } else {
490 #if 1
491 Bmat_addTo( A, iv, iw, em->NI[iv][r], em->NI[iw][s], em->A[iv][iw][r][s] );
492 #else
493 Bmat_addTo( A, iv, iw,
494 em->NI[iv][r], em->NI[iw][s],
495 em->A[iv][iw][r][s] );
496 #endif
497 }
498 } /* symmetries in the tangent form */
499 } /* assemble tangent matrix */
500
501 /* assemble mass matrix */
502 if (massKey > 0) {
503 /* handle symmetries in mass form */
504 if ( (sym==0) || ((sym==1)
505 &&(em->NI[iv][r]<=em->NI[iw][s])) ) {
506
507 /* one or both vr|vs are dirichlet */
508 if ( VDIRICHLET( em->TP[iv][r] )
509 || VDIRICHLET( em->TP[iw][s] ) ) {
510
511 /* vr is dirichlet */
512 if (VDIRICHLET( em->TP[iv][r] )) {
513 #if 1
514 Bmat_set( M,iv,iv, em->NI[iv][r], em->NI[iv][r], 1. );
515 #else
516 for (j=0; j<PDE_vec(thee->pde); j++)
517 Bmat_set( M,j,j,
518 em->NI[iv][r],
519 em->NI[iv][r], 1. );
520 #endif
521 }
522 /* vs is dirichlet */
523 if (VDIRICHLET( em->TP[iw][s] )) {
524 #if 1
525 Bmat_set( M,iw,iw, em->NI[iw][s], em->NI[iw][s], 1. );
526 #else
527 for (j=0; j<PDE_vec(thee->pde); j++)
528 Bmat_set( M,j,j,
529 em->NI[iw][s],
530 em->NI[iw][s], 1. );
531 #endif
532 }
533
534 /* neither vr|vs are dirichlet */
535 } else {
536 #if 1
537 Bmat_addTo( M, iv, iw, em->NI[iv][r], em->NI[iw][s], em->M[iv][iw][r][s] );
538 #else
539 Bmat_addTo( M, iv, iw,
540 em->NI[iv][r], em->NI[iw][s],
541 em->M[iv][iw][r][s] );
542 #endif
543 }
544 } /* symmetries in the mass form */
545 } /* assemble mass matrix */
546
547 } /* which W-basis function */
548 } /* s; loop over element matrix columns */
549 } /* if tangent or mass matrix needed */
550 } /* which V-basis function */
551 } /* r; loop over element matrix rows (and, faces of simplex) */
552 } /* if residual/load vector needed */
553 }
554
555 /*
556 * ***************************************************************************
557 * Routine: Aprx_assemEmat
558 *
559 * Purpose: Assemble one element matrix.
560 *
561 * Author: Michael Holst
562 * ***************************************************************************
563 */
Aprx_assemEmat(Aprx * thee,SS * sm,int evalKey,int energyKey,int residKey,int tangKey,int massKey,int bumpKey,Emat * em,Bvec * Wu,Bvec * Wud)564 VPUBLIC void Aprx_assemEmat(Aprx *thee, SS *sm,
565 int evalKey, int energyKey, int residKey, int tangKey, int massKey,
566 int bumpKey,
567 Emat *em,
568 Bvec *Wu, Bvec *Wud)
569 {
570 int face, someComponentIsNeumann, i, j;
571 TT t;
572
573 /* volume trans from master to this element (and back) */
574 Gem_buildVolumeTrans(thee->gm,sm,&t);
575
576 /* initialize the element matrix (zero it out) */
577 Aprx_initEmat(thee,bumpKey,&t,em);
578
579 /* once-per-element PDE initialization */
580 thee->pde->initElement(thee->pde,t.stype,t.gchart,t.vx,&t);
581
582 /* volume quadrature */
583 face = -1;
584 Aprx_quadEmat(thee,
585 evalKey,energyKey,residKey,tangKey,massKey,
586 bumpKey,
587 face,&t,em,Wu,Wud);
588
589 if(bumpKey) return;
590
591 /* now handle any neumann boundary faces this element may have */
592 for (face=0; face<Gem_dimVV(thee->gm); face++) {
593 t.focusFace = face;
594
595 someComponentIsNeumann = 0;
596 for (i=0; i<PDE_vec(thee->pde); i++) {
597 if (VNEUMANN( thee->pde->bmap[i][ SS_faceType(sm,face) ] )) {
598 someComponentIsNeumann = 1;
599 }
600 }
601
602 if ( someComponentIsNeumann ) {
603
604 /* boundary trans from master to this element (and back) */
605 Gem_buildSurfaceTrans(thee->gm,face,&t);
606
607 /* once-per-face PDE initialization */
608 thee->pde->initFace(thee->pde,SS_faceType(sm,face),
609 t.gchart,t.nvec[face]);
610
611 /* surface quadrature */
612 Aprx_quadEmat(thee,
613 evalKey,energyKey,residKey,tangKey,massKey,
614 bumpKey,
615 face,&t,em,Wu,Wud);
616 }
617 }
618 }
619
620 /*
621 * ***************************************************************************
622 * Routine: Aprx_assem
623 *
624 * Purpose: Assemble the tangent matrix and nonlinear residual
625 * vector (which reduces to the normal stiffness matrix and load
626 * vector in the case of a linear problem) for a general nonlinear
627 * system of m components in spatial dimension d.
628 *
629 * NOTES: The assembly modes for this routine are:
630 *
631 * evalKey == 0 ==> Primal problem: evaluation at z=[0+ud].
632 * == 1 ==> Primal problem: evaluation at z=[u+ud].
633 * == 2 ==> Dual problem: evaluation at z=[0+ud].
634 * == 3 ==> Dual problem: evaluation at z=[u+ud].
635 * == ? ==> Same as 0.
636 *
637 * energyKey == 0 ==> DON'T assemble an energy.
638 * == 1 ==> Assemble the energy J(z).
639 * == 2 ==> Assemble the energy 0.
640 * == ? ==> Same as 0.
641 *
642 * residKey == 0 ==> DON'T assemble a residual.
643 * == 1 ==> Assemble the residual F(z)(v).
644 * == 2 ==> Assemble the residual 0.
645 * == ? ==> Same as 0.
646 *
647 * tangKey == 0 ==> DON'T assemble a tangent matrix.
648 * == 1 ==> Assemble the tangent matrix DF(z)(w,v).
649 * == 2 ==> Assemble the tangent matrix 0.
650 * == ? ==> Same as 0.
651 *
652 * massKey == 0 ==> DON'T assemble a mass matrix.
653 * == 1 ==> Assemble the mass matrix p(w,v).
654 * == 2 ==> Assemble the mass matrix 0.
655 * == ? ==> Same as 0.
656 *
657 * NOTES: The nonlinear residual should always be ZERO at dirichlet nodes!
658 *
659 * Author: Michael Holst
660 * ***************************************************************************
661 */
Aprx_assem(Aprx * thee,int evalKey,int energyKey,int residKey,int tangKey,int massKey,int bumpKey,int ip[],double rp[],Bmat * A,Bmat * M,Bvec * F,Bvec * Wu,Bvec * Wud)662 VPUBLIC double Aprx_assem(Aprx *thee,
663 int evalKey, int energyKey, int residKey, int tangKey, int massKey,
664 int bumpKey,
665 int ip[], double rp[],
666 Bmat *A, Bmat *M, Bvec *F, Bvec *Wu, Bvec *Wud)
667 {
668 /* some variables */
669 int i, j, k, smid, vxid, vType, chart, cnt;
670 double energy, xm[3], Delt[MAXV];
671 SS *sm;
672 VV *vx;
673 Emat em;
674
675 /* do some i/o if this is a single assembly (e.g. linear problem) */
676 if ((evalKey==0) || (evalKey==2)) {
677 Vnm_tstart(30, "element assembly");
678 Vnm_print(0,"Aprx_assem: assembling..");
679 }
680
681 /* once-per-assembly PDE initialization */
682 thee->pde->initAssemble(thee->pde,ip,rp);
683
684 /* initialize the energy */
685 energy = 0.;
686
687 /* the BIG loop over all of the elements in the mesh.. */
688 smid = 0;
689 while ( (smid < Gem_numSS(thee->gm)) && (!Vnm_sigInt()) ) {
690 sm = Gem_SS(thee->gm,smid);
691
692 if ( (smid>0) && (smid % VPRTKEY) == 0 ) Vnm_print(0,"[AS:%d]",smid);
693
694 /* assemble the element matrix and the element vector */
695 Aprx_assemEmat(thee, sm,
696 evalKey, energyKey, residKey, tangKey, massKey,
697 bumpKey,
698 &em, Wu, Wud);
699
700 /* fan-out element matrix and element vector to global system */
701 Aprx_fanEmat(thee,
702 evalKey, energyKey, residKey, tangKey, massKey,
703 bumpKey,
704 &em, A, M, F);
705
706 /* accumulate the energy */
707 energy += em.J;
708
709 /* next sm */
710 smid++;
711
712 } /* loop over elements */
713
714 if(bumpKey) { /* bail out if we're solving the bump system */
715 /* do some i/o if this is a single assembly (linear problem) */
716 if ((evalKey==0) || (evalKey==2)) {
717 Vnm_print(0,"..done.\n");
718 Vnm_tstop(30, "assembly");
719 }
720
721 /* return the assembled energy */
722 return energy;
723 }
724
725 /* the BIG loop over all of the nodes in the mesh.. */
726 cnt = 0;
727 for (k=0; k<Bnode_numB(thee->B); k++) {
728 for (i=0; i<Bnode_numR(thee->B,k); i++) {
729 if ( (cnt>0) && (cnt % VPRTKEY) == 0 ) Vnm_print(0,"[AV:%d]",cnt);
730 cnt++;
731
732 /* MIKE: FIX THIS! */
733 vType = Bnode_data(thee->B,k)[i].type;
734 #if 1
735 chart = Bnode_data(thee->B,k)[i].chart;
736 #else
737 chart = VV_chart( Gem_VV(thee->gm,i) );
738 #endif
739 vxid = Bnode_data(thee->B,k)[i].myvert;
740
741 /* only deal with it if node is actually non-diri vertex of mesh */
742 if ( (!VDIRICHLET( vType )) && (vxid >= 0) && (vxid<Gem_numVV(thee->gm))) {
743 vx = Gem_VV(thee->gm,vxid);
744 for (j=0; j<Bnode_data(thee->B,k)[i].numx; j++)
745 xm[j] = Bnode_data(thee->B,k)[i].x[j];
746 for (j=0; j<PDE_vec(thee->pde); j++)
747 Delt[j] = 0.0;
748
749 /* deal with any delta-func contributions to the source term */
750 /* (assuming our basis functions are a Lagrange basis!) */
751 thee->pde->delta(thee->pde, vType, chart, xm, vx, Delt);
752
753 /* count ONLY the k-th component */
754 Bvec_addToB( F, k, vxid, Delt[k] );
755 }
756 }
757 }
758
759 /* build Bnode information into A and M */
760 Aprx_buildBRC(thee, A, M);
761
762 /* do some i/o if this is a single assembly (linear problem) */
763 if ((evalKey==0) || (evalKey==2)) {
764 Vnm_print(0,"..done.\n");
765 Vnm_tstop(30, "assembly");
766 }
767
768 /* return the assembled energy */
769 return energy;
770 }
771
772 /*
773 * ***************************************************************************
774 * Routine: Aprx_buildBRC
775 *
776 * Purpose: Build boundary Bnode information into two Bmats.
777 *
778 * Notes: Either/both Bmat objects can be VNULL without error;
779 * the result for that VNULL Bmat is a No-Op.
780 *
781 * Author: Michael Holst
782 * ***************************************************************************
783 */
Aprx_buildBRC(Aprx * thee,Bmat * A,Bmat * M)784 VPUBLIC void Aprx_buildBRC(Aprx *thee, Bmat *A, Bmat *M)
785 {
786 /* some variables */
787 int i, j;
788 int numBR[MAXV], numBC[MAXV], *BR[MAXV], *BC[MAXV];
789
790 /* sanity check: Bnode B better exist, since B/ibv created together */
791 VASSERT( thee->B != VNULL );
792
793 /* make boundary node arrays */
794 for (i=0; i<thee->numB; i++) {
795 numBR[i] = thee->numBV[i];
796 numBC[i] = thee->numBV[i];
797 if (numBR[i] > 0) {
798 BR[i] = Vmem_malloc(thee->vmem,numBR[i],sizeof(int));
799 }
800 if (numBC[i] > 0) {
801 BC[i] = Vmem_malloc(thee->vmem,numBC[i],sizeof(int));
802 }
803 for (j=0; j<numBR[i]; j++) {
804 BR[i][j] = thee->ibv[i][j];
805 }
806 for (j=0; j<numBC[i]; j++) {
807 BC[i][j] = thee->ibv[i][j];
808 }
809 }
810
811 /* copy boundary node arrays into A */
812 for (i=0; i<thee->numB; i++) {
813 for (j=0; j<thee->numB; j++) {
814 if (A != VNULL) {
815 Mat_buildBRC(A->AD[i][j], numBR[i], numBC[j], BR[i], BC[j]);
816 }
817 if (M != VNULL) {
818 Mat_buildBRC(M->AD[i][j], numBR[i], numBC[j], BR[i], BC[j]);
819 }
820 }
821 }
822
823 /* free boundary node arrays */
824 for (i=0; i<thee->numB; i++) {
825 if (numBR[i] > 0) {
826 Vmem_free(thee->vmem,numBR[i],sizeof(int),(void**)&BR[i]);
827 }
828 if (numBC[i] > 0) {
829 Vmem_free(thee->vmem,numBC[i],sizeof(int),(void**)&BC[i]);
830 }
831 }
832 }
833
834