APBS 3.0.0
Loading...
Searching...
No Matches
vpee.c
Go to the documentation of this file.
1
57#include "vpee.h"
58
59VPRIVATE int Vpee_userDefined(Vpee *thee,
60 SS *sm
61 );
62VPRIVATE int Vpee_ourSimp(Vpee *thee,
63 SS *sm,
64 int rcol
65 );
66VEXTERNC double Aprx_estNonlinResid(Aprx *thee,
67 SS *sm,
68 Bvec *u,
69 Bvec *ud,
70 Bvec *f
71 );
72VEXTERNC double Aprx_estLocalProblem(Aprx *thee,
73 SS *sm,
74 Bvec *u,
75 Bvec *ud,
76 Bvec *f);
77VEXTERNC double Aprx_estDualProblem(Aprx *thee,
78 SS *sm,
79 Bvec *u,
80 Bvec *ud,
81 Bvec *f
82 );
83
84/* ///////////////////////////////////////////////////////////////////////////
85// Class Vpee: Non-inlineable methods
87
88/* ///////////////////////////////////////////////////////////////////////////
89// Routine: Vpee_ctor
90//
91// Author: Nathan Baker
93VPUBLIC Vpee* Vpee_ctor(Gem *gm,
94 int localPartID,
95 int killFlag,
96 double killParam
97 ) {
98
99 Vpee *thee = VNULL;
100
101 /* Set up the structure */
102 thee = Vmem_malloc(VNULL, 1, sizeof(Vpee) );
103 VASSERT( thee != VNULL);
104 VASSERT( Vpee_ctor2(thee, gm, localPartID, killFlag, killParam));
105
106 return thee;
107}
108
109/* ///////////////////////////////////////////////////////////////////////////
110// Routine: Vpee_ctor2
111//
112// Author: Nathan Baker
114VPUBLIC int Vpee_ctor2(Vpee *thee,
115 Gem *gm,
116 int localPartID,
117 int killFlag,
118 double killParam
119 ) {
120
121 int ivert,
122 nLocalVerts;
123 SS *simp;
124 VV *vert;
125 double radius,
126 dx,
127 dy,
128 dz;
129
130 VASSERT(thee != VNULL);
131
132 /* Sanity check on input values */
133 if (killFlag == 0) {
134 Vnm_print(0, "Vpee_ctor2: No error attenuation outside partition.\n");
135 } else if (killFlag == 1) {
136 Vnm_print(0, "Vpee_ctor2: Error outside local partition ignored.\n");
137 } else if (killFlag == 2) {
138 Vnm_print(0, "Vpee_ctor2: Error ignored outside sphere with radius %4.3f times the radius of the circumscribing sphere\n", killParam);
139 if (killParam < 1.0) {
140 Vnm_print(2, "Vpee_ctor2: Warning! Parameter killParam = %4.3 < 1.0!\n",
141 killParam);
142 Vnm_print(2, "Vpee_ctor2: This may result in non-optimal marking and refinement!\n");
143 }
144 } else if (killFlag == 3) {
145 Vnm_print(0, "Vpee_ctor2: Error outside local partition and immediate neighbors ignored [NOT IMPLEMENTED].\n");
146 } else {
147 Vnm_print(2, "Vpee_ctor2: UNRECOGNIZED killFlag PARAMETER! BAILING!.\n");
148 VASSERT(0);
149 }
150
151 thee->gm = gm;
152 thee->localPartID = localPartID;
153 thee->killFlag = killFlag;
154 thee->killParam = killParam;
155 thee->mem = Vmem_ctor("APBS::VPEE");
156
157 /* Now, figure out the center of geometry for the local partition. The
158 * general plan is to loop through the vertices, loop through the
159 * vertices' simplex lists and find the vertices with simplices containing
160 * chart values we're interested in. */
161 thee->localPartCenter[0] = 0.0;
162 thee->localPartCenter[1] = 0.0;
163 thee->localPartCenter[2] = 0.0;
164 nLocalVerts = 0;
165 for (ivert=0; ivert<Gem_numVV(thee->gm); ivert++) {
166 vert = Gem_VV(thee->gm, ivert);
167 simp = VV_firstSS(vert);
168 VASSERT(simp != VNULL);
169 while (simp != VNULL) {
170 if (SS_chart(simp) == thee->localPartID) {
171 thee->localPartCenter[0] += VV_coord(vert, 0);
172 thee->localPartCenter[1] += VV_coord(vert, 1);
173 thee->localPartCenter[2] += VV_coord(vert, 2);
174 nLocalVerts++;
175 break;
176 }
177 simp = SS_link(simp, vert);
178 }
179 }
180 VASSERT(nLocalVerts > 0);
181 thee->localPartCenter[0] =
182 thee->localPartCenter[0]/((double)(nLocalVerts));
183 thee->localPartCenter[1] =
184 thee->localPartCenter[1]/((double)(nLocalVerts));
185 thee->localPartCenter[2] =
186 thee->localPartCenter[2]/((double)(nLocalVerts));
187 Vnm_print(0, "Vpee_ctor2: Part %d centered at (%4.3f, %4.3f, %4.3f)\n",
188 thee->localPartID, thee->localPartCenter[0], thee->localPartCenter[1],
189 thee->localPartCenter[2]);
190
191
192 /* Now, figure out the radius of the sphere circumscribing the local
193 * partition. We need to keep track of vertices so we don't double count
194 * them. */
195 thee->localPartRadius = 0.0;
196 for (ivert=0; ivert<Gem_numVV(thee->gm); ivert++) {
197 vert = Gem_VV(thee->gm, ivert);
198 simp = VV_firstSS(vert);
199 VASSERT(simp != VNULL);
200 while (simp != VNULL) {
201 if (SS_chart(simp) == thee->localPartID) {
202 dx = thee->localPartCenter[0] - VV_coord(vert, 0);
203 dy = thee->localPartCenter[1] - VV_coord(vert, 1);
204 dz = thee->localPartCenter[2] - VV_coord(vert, 2);
205 radius = dx*dx + dy*dy + dz*dz;
206 if (radius > thee->localPartRadius) thee->localPartRadius =
207 radius;
208 break;
209 }
210 simp = SS_link(simp, vert);
211 }
212 }
213 thee->localPartRadius = VSQRT(thee->localPartRadius);
214 Vnm_print(0, "Vpee_ctor2: Part %d has circumscribing sphere of radius %4.3f\n",
215 thee->localPartID, thee->localPartRadius);
216
217 return 1;
218}
219
220/* ///////////////////////////////////////////////////////////////////////////
221// Routine: Vpee_dtor
222//
223// Author: Nathan Baker
225VPUBLIC void Vpee_dtor(Vpee **thee) {
226
227 if ((*thee) != VNULL) {
228 Vpee_dtor2(*thee);
229 Vmem_free(VNULL, 1, sizeof(Vpee), (void **)thee);
230 (*thee) = VNULL;
231 }
232
233}
234
235/* ///////////////////////////////////////////////////////////////////////////
236// Routine: Vpee_dtor2
237//
238// Author: Nathan Baker
240VPUBLIC void Vpee_dtor2(Vpee *thee) {
241 Vmem_dtor(&(thee->mem));
242}
243
244/* ///////////////////////////////////////////////////////////////////////////
245// Routine: Vpee_markRefine
246//
247// Author: Nathan Baker (and Michael Holst: the author of AM_markRefine, on
248// which this is based)
250VPUBLIC int Vpee_markRefine(Vpee *thee,
251 AM *am,
252 int level,
253 int akey,
254 int rcol,
255 double etol,
256 int bkey
257 ) {
258
259 Aprx *aprx;
260 int marked = 0,
261 markMe,
262 i,
263 smid,
264 count,
265 currentQ;
266 double minError = 0.0,
267 maxError = 0.0,
268 errEst = 0.0,
269 mlevel,
270 barrier;
271 SS *sm;
272
273
274 VASSERT(thee != VNULL);
275
276 /* Get the Aprx object from AM */
277 aprx = am->aprx;
278
279 /* input check and some i/o */
280 if ( ! ((-1 <= akey) && (akey <= 4)) ) {
281 Vnm_print(0,"Vpee_markRefine: bad refine key; simplices marked = %d\n",
282 marked);
283 return marked;
284 }
285
286 /* For uniform markings, we have no effect */
287 if ((-1 <= akey) && (akey <= 0)) {
288 marked = Gem_markRefine(thee->gm, akey, rcol);
289 return marked;
290 }
291
292 /* Informative I/O */
293 if (akey == 2) {
294 Vnm_print(0,"Vpee_estRefine: using Aprx_estNonlinResid().\n");
295 } else if (akey == 3) {
296 Vnm_print(0,"Vpee_estRefine: using Aprx_estLocalProblem().\n");
297 } else if (akey == 4) {
298 Vnm_print(0,"Vpee_estRefine: using Aprx_estDualProblem().\n");
299 } else {
300 Vnm_print(0,"Vpee_estRefine: bad key given; simplices marked = %d\n",
301 marked);
302 return marked;
303 }
304 if (thee->killFlag == 0) {
305 Vnm_print(0, "Vpee_markRefine: No error attenuation -- simplices in all partitions will be marked.\n");
306 } else if (thee->killFlag == 1) {
307 Vnm_print(0, "Vpee_markRefine: Maximum error attenuation -- only simplices in local partition will be marked.\n");
308 } else if (thee->killFlag == 2) {
309 Vnm_print(0, "Vpee_markRefine: Spherical error attenutation -- simplices within a sphere of %4.3f times the size of the partition will be marked\n",
310 thee->killParam);
311 } else if (thee->killFlag == 2) {
312 Vnm_print(0, "Vpee_markRefine: Neighbor-based error attenuation -- simplices in the local and neighboring partitions will be marked [NOT IMPLEMENTED]!\n");
313 VASSERT(0);
314 } else {
315 Vnm_print(2,"Vpee_markRefine: bogus killFlag given; simplices marked = %d\n",
316 marked);
317 return marked;
318 }
319
320 /* set the barrier type */
321 mlevel = (etol*etol) / Gem_numSS(thee->gm);
322 if (bkey == 0) {
323 barrier = (etol*etol);
324 Vnm_print(0,"Vpee_estRefine: forcing [err per S] < [TOL] = %g\n",
325 barrier);
326 } else if (bkey == 1) {
327 barrier = mlevel;
328 Vnm_print(0,"Vpee_estRefine: forcing [err per S] < [(TOL^2/numS)^{1/2}] = %g\n",
329 VSQRT(barrier));
330 } else {
331 Vnm_print(0,"Vpee_estRefine: bad bkey given; simplices marked = %d\n",
332 marked);
333 return marked;
334 }
335
336 /* timer */
337 Vnm_tstart(30, "error estimation");
338
339 /* count = num generations to produce from marked simplices (minimally) */
340 count = 1; /* must be >= 1 */
341
342 /* check the refinement Q for emptyness */
343 currentQ = 0;
344 if (Gem_numSQ(thee->gm,currentQ) > 0) {
345 Vnm_print(0,"Vpee_markRefine: non-empty refinement Q%d....clearing..",
346 currentQ);
347 Gem_resetSQ(thee->gm,currentQ);
348 Vnm_print(0,"..done.\n");
349 }
350 if (Gem_numSQ(thee->gm,!currentQ) > 0) {
351 Vnm_print(0,"Vpee_markRefine: non-empty refinement Q%d....clearing..",
352 !currentQ);
353 Gem_resetSQ(thee->gm,!currentQ);
354 Vnm_print(0,"..done.\n");
355 }
356 VASSERT( Gem_numSQ(thee->gm,currentQ) == 0 );
357 VASSERT( Gem_numSQ(thee->gm,!currentQ) == 0 );
358
359 /* clear everyone's refinement flags */
360 Vnm_print(0,"Vpee_markRefine: clearing all simplex refinement flags..");
361 for (i=0; i<Gem_numSS(thee->gm); i++) {
362 if ( (i>0) && (i % VPRTKEY) == 0 ) Vnm_print(0,"[MS:%d]",i);
363 sm = Gem_SS(thee->gm,i);
364 SS_setRefineKey(sm,currentQ,0);
365 SS_setRefineKey(sm,!currentQ,0);
366 SS_setRefinementCount(sm,0);
367 }
368 Vnm_print(0,"..done.\n");
369
370 /* NON-ERROR-BASED METHODS */
371 /* Simplex flag clearing */
372 if (akey == -1) return marked;
373 /* Uniform & user-defined refinement*/
374 if ((akey == 0) || (akey == 1)) {
375 smid = 0;
376 while ( smid < Gem_numSS(thee->gm)) {
377 /* Get the simplex and find out if it's markable */
378 sm = Gem_SS(thee->gm,smid);
379 markMe = Vpee_ourSimp(thee, sm, rcol);
380 if (markMe) {
381 if (akey == 0) {
382 marked++;
383 Gem_appendSQ(thee->gm,currentQ, sm);
384 SS_setRefineKey(sm,currentQ,1);
385 SS_setRefinementCount(sm,count);
386 } else if (Vpee_userDefined(thee, sm)) {
387 marked++;
388 Gem_appendSQ(thee->gm,currentQ, sm);
389 SS_setRefineKey(sm,currentQ,1);
390 SS_setRefinementCount(sm,count);
391 }
392 }
393 smid++;
394 }
395 }
396
397 /* ERROR-BASED METHODS */
398 /* gerror = global error accumulation */
399 aprx->gerror = 0.;
400
401 /* traverse the simplices and process the error estimates */
402 Vnm_print(0,"Vpee_markRefine: estimating error..");
403 smid = 0;
404 while ( smid < Gem_numSS(thee->gm)) {
405
406 /* Get the simplex and find out if it's markable */
407 sm = Gem_SS(thee->gm,smid);
408 markMe = Vpee_ourSimp(thee, sm, rcol);
409
410 if ( (smid>0) && (smid % VPRTKEY) == 0 ) Vnm_print(0,"[MS:%d]",smid);
411
412 /* Produce an error estimate for this element if it is in the set */
413 if (markMe) {
414 if (akey == 2) {
415 errEst = Aprx_estNonlinResid(aprx, sm, am->u,am->ud,am->f);
416 } else if (akey == 3) {
417 errEst = Aprx_estLocalProblem(aprx, sm, am->u,am->ud,am->f);
418 } else if (akey == 4) {
419 errEst = Aprx_estDualProblem(aprx, sm, am->u,am->ud,am->f);
420 }
421 VASSERT( errEst >= 0. );
422
423 /* if error estimate above tol, mark element for refinement */
424 if ( errEst > barrier ) {
425 marked++;
426 Gem_appendSQ(thee->gm,currentQ, sm); /*add to refinement Q*/
427 SS_setRefineKey(sm,currentQ,1); /* note now on refine Q */
428 SS_setRefinementCount(sm,count); /* refine X many times? */
429 }
430
431 /* keep track of min/max errors over the mesh */
432 minError = VMIN2( VSQRT(VABS(errEst)), minError );
433 maxError = VMAX2( VSQRT(VABS(errEst)), maxError );
434
435 /* store the estimate */
436 Bvec_set( aprx->wev, smid, errEst );
437
438 /* accumlate into global error (errEst is SQUAREd already) */
439 aprx->gerror += errEst;
440
441 /* otherwise store a zero for the estimate */
442 } else {
443 Bvec_set( aprx->wev, smid, 0. );
444 }
445
446 smid++;
447 }
448
449 /* do some i/o */
450 Vnm_print(0,"..done. [marked=<%d/%d>]\n",marked,Gem_numSS(thee->gm));
451 Vnm_print(0,"Vpee_estRefine: TOL=<%g> Global_Error=<%g>\n",
452 etol, aprx->gerror);
453 Vnm_print(0,"Vpee_estRefine: (TOL^2/numS)^{1/2}=<%g> Max_Ele_Error=<%g>\n",
454 VSQRT(mlevel),maxError);
455 Vnm_tstop(30, "error estimation");
456
457 /* check for making the error tolerance */
458 if ((bkey == 1) && (aprx->gerror <= etol)) {
459 Vnm_print(0,
460 "Vpee_estRefine: *********************************************\n");
461 Vnm_print(0,
462 "Vpee_estRefine: Global Error criterion met; setting marked=0.\n");
463 Vnm_print(0,
464 "Vpee_estRefine: *********************************************\n");
465 marked = 0;
466 }
467
468
469 /* return */
470 return marked;
471
472}
473
474/* ///////////////////////////////////////////////////////////////////////////
475// Routine: Vpee_numSS
476//
477// Author: Nathan Baker
479VPUBLIC int Vpee_numSS(Vpee *thee) {
480 int num = 0;
481 int isimp;
482
483 for (isimp=0; isimp<Gem_numSS(thee->gm); isimp++) {
484 if (SS_chart(Gem_SS(thee->gm, isimp)) == thee->localPartID) num++;
485 }
486
487 return num;
488}
489
490/* ///////////////////////////////////////////////////////////////////////////
491// Routine: Vpee_userDefined
492//
493// Purpose: Reduce code bloat by wrapping up the common steps for getting the
494// user-defined error estimate
495//
496// Author: Nathan Baker
498VPRIVATE int Vpee_userDefined(Vpee *thee,
499 SS *sm
500 ) {
501
502 int ivert,
503 icoord,
504 chart[4],
505 fType[4],
506 vType[4];
507 double vx[4][3];
508
509 for (ivert=0; ivert<Gem_dimVV(thee->gm); ivert++) {
510 fType[ivert] = SS_faceType(sm,ivert);
511 vType[ivert] = VV_type(SS_vertex(sm,ivert) );
512 chart[ivert] = VV_chart(SS_vertex(sm,ivert) );
513 for (icoord=0; icoord<Gem_dimII(thee->gm); icoord++) {
514 vx[ivert][icoord] = VV_coord(SS_vertex(sm,ivert), icoord );
515 }
516 }
517 return thee->gm->pde->markSimplex(Gem_dim(thee->gm), Gem_dimII(thee->gm),
518 SS_type(sm), fType, vType, chart, vx, sm);
519}
520
521/* ///////////////////////////////////////////////////////////////////////////
522// Routine: Vpee_ourSimp
523//
524// Purpose: Reduce code bloat by wrapping up the common steps for determining
525// whether the given simplex can be marked (i.e., belongs to our
526// partition or overlap region)
527//
528// Returns: 1 if could be marked, 0 otherwise
529//
530// Author: Nathan Baker
532VPRIVATE int Vpee_ourSimp(Vpee *thee,
533 SS *sm,
534 int rcol
535 ) {
536
537 int ivert;
538 double dist,
539 dx,
540 dy,
541 dz;
542
543 if (thee->killFlag == 0) return 1;
544 else if (thee->killFlag == 1) {
545 if ((SS_chart(sm) == rcol) || (rcol < 0)) return 1;
546 } else if (thee->killFlag == 2) {
547 if (rcol < 0) return 1;
548 else {
549 /* We can only do distance-based searches on the local partition */
550 VASSERT(rcol == thee->localPartID);
551 /* Find the closest distance between this simplex and the
552 * center of the local partition and check it against
553 * (thee->localPartRadius*thee->killParam) */
554 dist = 0;
555 for (ivert=0; ivert<SS_dimVV(sm); ivert++) {
556 dx = VV_coord(SS_vertex(sm, ivert), 0) -
557 thee->localPartCenter[0];
558 dy = VV_coord(SS_vertex(sm, ivert), 1) -
559 thee->localPartCenter[1];
560 dz = VV_coord(SS_vertex(sm, ivert), 2) -
561 thee->localPartCenter[2];
562 dist = VSQRT((dx*dx + dy*dy + dz*dz));
563 }
564 if (dist < thee->localPartRadius*thee->killParam) return 1;
565 }
566 } else if (thee->killFlag == 3) VASSERT(0);
567 else VASSERT(0);
568
569 return 0;
570
571}
VPUBLIC int Vpee_ctor2(Vpee *thee, Gem *gm, int localPartID, int killFlag, double killParam)
FORTRAN stub to construct the Vpee object.
Definition vpee.c:114
Contains public data members for Vpee class/module.
Definition vpee.h:89
Contains declarations for class Vpee.