APBS  1.5
vpee.c
Go to the documentation of this file.
1 
57 #include "vpee.h"
58 
59 VPRIVATE int Vpee_userDefined(Vpee *thee,
60  SS *sm
61  );
62 VPRIVATE int Vpee_ourSimp(Vpee *thee,
63  SS *sm,
64  int rcol
65  );
66 VEXTERNC double Aprx_estNonlinResid(Aprx *thee,
67  SS *sm,
68  Bvec *u,
69  Bvec *ud,
70  Bvec *f
71  );
72 VEXTERNC double Aprx_estLocalProblem(Aprx *thee,
73  SS *sm,
74  Bvec *u,
75  Bvec *ud,
76  Bvec *f);
77 VEXTERNC 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
93 VPUBLIC 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
114 VPUBLIC 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
225 VPUBLIC 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
240 VPUBLIC 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)
250 VPUBLIC 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
479 VPUBLIC 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
498 VPRIVATE 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
532 VPRIVATE 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 }
Contains declarations for class Vpee.
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