Actual source code: rgring.c

slepc-3.11.2 2019-07-30
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2019, Universitat Politecnica de Valencia, Spain

  6:    This file is part of SLEPc.
  7:    SLEPc is distributed under a 2-clause BSD license (see LICENSE).
  8:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  9: */
 10: /*
 11:    Ring region, similar to the ellipse but with a start and end angle,
 12:    together with the width
 13: */

 15: #include <slepc/private/rgimpl.h>      /*I "slepcrg.h" I*/
 16: #include <petscdraw.h>

 18: typedef struct {
 19:   PetscScalar center;     /* center of the ellipse */
 20:   PetscReal   radius;     /* radius of the ellipse */
 21:   PetscReal   vscale;     /* vertical scale of the ellipse */
 22:   PetscReal   start_ang;  /* start angle */
 23:   PetscReal   end_ang;    /* end angle */
 24:   PetscReal   width;      /* ring width */
 25: } RG_RING;

 27: static PetscErrorCode RGRingSetParameters_Ring(RG rg,PetscScalar center,PetscReal radius,PetscReal vscale,PetscReal start_ang,PetscReal end_ang,PetscReal width)
 28: {
 29:   RG_RING *ctx = (RG_RING*)rg->data;

 32:   ctx->center = center;
 33:   if (radius == PETSC_DEFAULT) {
 34:     ctx->radius = 1.0;
 35:   } else {
 36:     if (radius<=0.0) SETERRQ(PetscObjectComm((PetscObject)rg),PETSC_ERR_ARG_OUTOFRANGE,"The radius argument must be > 0.0");
 37:     ctx->radius = radius;
 38:   }
 39:   if (vscale<=0.0) SETERRQ(PetscObjectComm((PetscObject)rg),PETSC_ERR_ARG_OUTOFRANGE,"The vscale argument must be > 0.0");
 40:   ctx->vscale = vscale;
 41:   if (start_ang == PETSC_DEFAULT) {
 42:     ctx->start_ang = 0.0;
 43:   } else {
 44:     if (start_ang<0.0) SETERRQ(PetscObjectComm((PetscObject)rg),PETSC_ERR_ARG_OUTOFRANGE,"The right-hand side angle argument must be >= 0.0");
 45:     if (start_ang>1.0) SETERRQ(PetscObjectComm((PetscObject)rg),PETSC_ERR_ARG_OUTOFRANGE,"The right-hand side angle argument must be <= 1.0");
 46:     ctx->start_ang = start_ang;
 47:   }
 48:   if (end_ang == PETSC_DEFAULT) {
 49:     ctx->end_ang = 1.0;
 50:   } else {
 51:     if (end_ang<0.0) SETERRQ(PetscObjectComm((PetscObject)rg),PETSC_ERR_ARG_OUTOFRANGE,"The left-hand side angle argument must be >= 0.0");
 52:     if (end_ang>1.0) SETERRQ(PetscObjectComm((PetscObject)rg),PETSC_ERR_ARG_OUTOFRANGE,"The left-hand side angle argument must be <= 1.0");
 53:     ctx->end_ang = end_ang;
 54:   }
 55: #if !defined(PETSC_USE_COMPLEX)
 56:   if (ctx->start_ang+ctx->end_ang!=1.0) SETERRQ(PetscObjectComm((PetscObject)rg),PETSC_ERR_ARG_WRONG,"In real scalars the region must be symmetric wrt real axis");
 57: #endif
 58:   if (width == PETSC_DEFAULT) {
 59:     ctx->width = 0.1;
 60:   } else {
 61:     if (width<=0.0) SETERRQ(PetscObjectComm((PetscObject)rg),PETSC_ERR_ARG_OUTOFRANGE,"The width argument must be > 0.0");
 62:     ctx->width = width;
 63:   }
 64:   return(0);
 65: }

 67: /*@
 68:    RGRingSetParameters - Sets the parameters defining the ring region.

 70:    Logically Collective on RG

 72:    Input Parameters:
 73: +  rg        - the region context
 74: .  center    - center of the ellipse
 75: .  radius    - radius of the ellipse
 76: .  vscale    - vertical scale of the ellipse
 77: .  start_ang - the right-hand side angle
 78: .  end_ang   - the left-hand side angle
 79: -  width     - width of the ring

 81:    Options Database Keys:
 82: +  -rg_ring_center     - Sets the center
 83: .  -rg_ring_radius     - Sets the radius
 84: .  -rg_ring_vscale     - Sets the vertical scale
 85: .  -rg_ring_startangle - Sets the right-hand side angle
 86: .  -rg_ring_endangle   - Sets the left-hand side angle
 87: -  -rg_ring_width      - Sets the width of the ring

 89:    Notes:
 90:    The values of center, radius and vscale have the same meaning as in the
 91:    ellipse region. The startangle and endangle define the span of the ring
 92:    (by default it is the whole ring), while the width is the separation
 93:    between the two concentric ellipses (above and below the radius by
 94:    width/2).

 96:    The start and end angles are expressed as a fraction of the circumference:
 97:    the allowed range is [0..1], with 0 corresponding to 0 radians, 0.25 to
 98:    pi/2 radians, and so on. It is allowed to have startangle>endangle, in
 99:    which case the ring region crosses over the zero angle.

101:    In the case of complex scalars, a complex center can be provided in the
102:    command line with [+/-][realnumber][+/-]realnumberi with no spaces, e.g.
103:    -rg_ring_center 1.0+2.0i

105:    When PETSc is built with real scalars, the center is restricted to a real value,
106:    and the start and end angles must be such that the region is symmetric with
107:    respect to the real axis.

109:    Level: advanced

111: .seealso: RGRingGetParameters()
112: @*/
113: PetscErrorCode RGRingSetParameters(RG rg,PetscScalar center,PetscReal radius,PetscReal vscale,PetscReal start_ang,PetscReal end_ang,PetscReal width)
114: {

125:   PetscTryMethod(rg,"RGRingSetParameters_C",(RG,PetscScalar,PetscReal,PetscReal,PetscReal,PetscReal,PetscReal),(rg,center,radius,vscale,start_ang,end_ang,width));
126:   return(0);
127: }

129: static PetscErrorCode RGRingGetParameters_Ring(RG rg,PetscScalar *center,PetscReal *radius,PetscReal *vscale,PetscReal *start_ang,PetscReal *end_ang,PetscReal *width)
130: {
131:   RG_RING *ctx = (RG_RING*)rg->data;

134:   if (center)    *center    = ctx->center;
135:   if (radius)    *radius    = ctx->radius;
136:   if (vscale)    *vscale    = ctx->vscale;
137:   if (start_ang) *start_ang = ctx->start_ang;
138:   if (end_ang)   *end_ang   = ctx->end_ang;
139:   if (width)     *width     = ctx->width;
140:   return(0);
141: }

143: /*@
144:    RGRingGetParameters - Gets the parameters that define the ring region.

146:    Not Collective

148:    Input Parameter:
149: .  rg     - the region context

151:    Output Parameters:
152: +  center    - center of the region
153: .  radius    - radius of the region
154: .  vscale    - vertical scale of the region
155: .  start_ang - the right-hand side angle
156: .  end_ang   - the left-hand side angle
157: -  width     - width of the ring

159:    Level: advanced

161: .seealso: RGRingSetParameters()
162: @*/
163: PetscErrorCode RGRingGetParameters(RG rg,PetscScalar *center,PetscReal *radius,PetscReal *vscale,PetscReal *start_ang,PetscReal *end_ang,PetscReal *width)
164: {

169:   PetscUseMethod(rg,"RGRingGetParameters_C",(RG,PetscScalar*,PetscReal*,PetscReal*,PetscReal*,PetscReal*,PetscReal*),(rg,center,radius,vscale,start_ang,end_ang,width));
170:   return(0);
171: }

173: PetscErrorCode RGView_Ring(RG rg,PetscViewer viewer)
174: {
176:   RG_RING        *ctx = (RG_RING*)rg->data;
177:   int            winw,winh;
178:   PetscBool      isdraw,isascii;
179:   PetscDraw      draw;
180:   PetscDrawAxis  axis;
181:   PetscReal      cx,cy,radius,width,ab,cd,lx,ly,w,end_ang,x1,y1,x2,y2,r,theta,scale=1.2;
182:   char           str[50];

185:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
186:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
187:   if (isascii) {
188:     SlepcSNPrintfScalar(str,50,ctx->center,PETSC_FALSE);
189:     PetscViewerASCIIPrintf(viewer,"  center: %s, radius: %g, vscale: %g, start angle: %g, end angle: %g, ring width: %g\n",str,RGShowReal(ctx->radius),RGShowReal(ctx->vscale),(double)ctx->start_ang,(double)ctx->end_ang,(double)ctx->width);
190:   } else if (isdraw) {
191:     PetscViewerDrawGetDraw(viewer,0,&draw);
192:     PetscDrawCheckResizedWindow(draw);
193:     PetscDrawGetWindowSize(draw,&winw,&winh);
194:     winw = PetscMax(winw,1); winh = PetscMax(winh,1);
195:     PetscDrawClear(draw);
196:     PetscDrawSetTitle(draw,"Ring region");
197:     PetscDrawAxisCreate(draw,&axis);
198:     cx = PetscRealPart(ctx->center)*rg->sfactor;
199:     cy = PetscImaginaryPart(ctx->center)*rg->sfactor;
200:     radius = ctx->radius*rg->sfactor;
201:     width  = ctx->width*rg->sfactor;
202:     lx = 2*(radius+width);
203:     ly = 2*(radius+width)*ctx->vscale;
204:     ab = cx;
205:     cd = cy;
206:     w  = scale*PetscMax(lx/winw,ly/winh)/2;
207:     PetscDrawAxisSetLimits(axis,ab-w*winw,ab+w*winw,cd-w*winh,cd+w*winh);
208:     PetscDrawAxisDraw(axis);
209:     PetscDrawAxisDestroy(&axis);
210:     /* draw outer ellipse */
211:     PetscDrawEllipse(draw,cx,cy,2*(radius+width),2*(radius+width)*ctx->vscale,PETSC_DRAW_ORANGE);
212:     /* remove inner part */
213:     PetscDrawEllipse(draw,cx,cy,2*(radius-width),2*(radius-width)*ctx->vscale,PETSC_DRAW_WHITE);
214:     if (ctx->start_ang!=ctx->end_ang) {
215:       /* remove section from end_ang to start_ang */
216:       end_ang = (ctx->start_ang<ctx->end_ang)? ctx->end_ang-1: ctx->end_ang;
217:       theta = end_ang;
218:       r = scale*(radius+width);
219:       if (ctx->vscale>1) r *= ctx->vscale;
220:       x1 = PetscMin(PetscMax(ab+r*PetscCosReal(2.0*PETSC_PI*theta),ab-w*winw),ab+w*winw);
221:       y1 = PetscMin(PetscMax(cd+r*PetscSinReal(2.0*PETSC_PI*theta),cd-w*winh),cd+w*winh);
222:       do {
223:         theta = PetscMin(PetscFloorReal(8*theta+1)/8,ctx->start_ang);
224:         x2 = PetscMin(PetscMax(ab+r*PetscCosReal(2.0*PETSC_PI*theta),ab-w*winw),ab+w*winw);
225:         y2 = PetscMin(PetscMax(cd+r*PetscSinReal(2.0*PETSC_PI*theta),cd-w*winh),cd+w*winh);
226:         PetscDrawTriangle(draw,cx,cy,x1,y1,x2,y2,PETSC_DRAW_WHITE,PETSC_DRAW_WHITE,PETSC_DRAW_WHITE);
227:         x1 = x2; y1 = y2;
228:       } while (theta<ctx->start_ang);
229:     }
230:     PetscDrawFlush(draw);
231:     PetscDrawSave(draw);
232:     PetscDrawPause(draw);
233:   }
234:   return(0);
235: }

237: PetscErrorCode RGIsTrivial_Ring(RG rg,PetscBool *trivial)
238: {
239:   RG_RING *ctx = (RG_RING*)rg->data;

242:   if (rg->complement) *trivial = PetscNot(ctx->radius);
243:   else *trivial = PetscNot(ctx->radius<PETSC_MAX_REAL);
244:   return(0);
245: }

247: PetscErrorCode RGComputeContour_Ring(RG rg,PetscInt n,PetscScalar *cr,PetscScalar *ci)
248: {
249:   RG_RING   *ctx = (RG_RING*)rg->data;
250:   PetscReal theta,start_ang;
251:   PetscInt  i,n2=n/2;

254:   start_ang = (ctx->start_ang>ctx->end_ang)? ctx->start_ang-1: ctx->start_ang;
255:   for (i=0;i<n;i++) {
256:     if (i < n2) {
257:       theta = ((ctx->end_ang-start_ang)*i/n2 + start_ang)*2.0*PETSC_PI;
258: #if defined(PETSC_USE_COMPLEX)
259:       cr[i] = ctx->center + (ctx->radius+ctx->width/2.0)*(PetscCosReal(theta)+ctx->vscale*PetscSinReal(theta)*PETSC_i);
260: #else
261:       cr[i] = ctx->center + (ctx->radius+ctx->width/2.0)*PetscCosReal(theta);
262:       ci[i] = (ctx->radius+ctx->width/2.0)*ctx->vscale*PetscSinReal(theta);
263: #endif
264:     } else {
265:       theta = ((ctx->end_ang-start_ang)*(n-i)/n2 + start_ang)*2.0*PETSC_PI;
266: #if defined(PETSC_USE_COMPLEX)
267:       cr[i] = ctx->center + (ctx->radius-ctx->width/2.0)*(PetscCosReal(theta)+ctx->vscale*PetscSinReal(theta)*PETSC_i);
268: #else
269:       cr[i] = ctx->center + (ctx->radius-ctx->width/2.0)*PetscCosReal(theta);
270:       ci[i] = (ctx->radius-ctx->width/2.0)*ctx->vscale*PetscSinReal(theta);
271: #endif
272:     }
273:   }
274:   return(0);
275: }

277: PetscErrorCode RGComputeBoundingBox_Ring(RG rg,PetscReal *a,PetscReal *b,PetscReal *c,PetscReal *d)
278: {
279:   RG_RING *ctx = (RG_RING*)rg->data;

282:   /* current implementation does not return a tight bounding box */
283:   if (a) *a = PetscRealPart(ctx->center) - (ctx->radius+ctx->width/2.0);
284:   if (b) *b = PetscRealPart(ctx->center) + (ctx->radius+ctx->width/2.0);
285:   if (c) *c = PetscImaginaryPart(ctx->center) - (ctx->radius+ctx->width/2.0)*ctx->vscale;
286:   if (d) *d = PetscImaginaryPart(ctx->center) + (ctx->radius+ctx->width/2.0)*ctx->vscale;
287:   return(0);
288: }

290: PetscErrorCode RGCheckInside_Ring(RG rg,PetscReal px,PetscReal py,PetscInt *inside)
291: {
292:   RG_RING   *ctx = (RG_RING*)rg->data;
293:   PetscReal dx,dy,r;

296:   /* outer ellipse */
297: #if defined(PETSC_USE_COMPLEX)
298:   dx = (px-PetscRealPart(ctx->center))/(ctx->radius+ctx->width/2.0);
299:   dy = (py-PetscImaginaryPart(ctx->center))/(ctx->radius+ctx->width/2.0);
300: #else
301:   dx = (px-ctx->center)/(ctx->radius+ctx->width/2.0);
302:   dy = py/(ctx->radius+ctx->width/2.0);
303: #endif
304:   r = 1.0-dx*dx-(dy*dy)/(ctx->vscale*ctx->vscale);
305:   *inside = PetscSign(r);
306:   /* inner ellipse */
307: #if defined(PETSC_USE_COMPLEX)
308:   dx = (px-PetscRealPart(ctx->center))/(ctx->radius-ctx->width/2.0);
309:   dy = (py-PetscImaginaryPart(ctx->center))/(ctx->radius-ctx->width/2.0);
310: #else
311:   dx = (px-ctx->center)/(ctx->radius-ctx->width/2.0);
312:   dy = py/(ctx->radius-ctx->width/2.0);
313: #endif
314:   r = -1.0+dx*dx+(dy*dy)/(ctx->vscale*ctx->vscale);
315:   *inside *= PetscSign(r);
316:   if (*inside == 1) {  /* check angles */
317: #if defined(PETSC_USE_COMPLEX)
318:     dx = (px-PetscRealPart(ctx->center));
319:     dy = (py-PetscImaginaryPart(ctx->center));
320: #else
321:     dx = px-ctx->center;
322:     dy = py;
323: #endif
324:     if (dx == 0) {
325:       if (dy == 0) r = -1;
326:       else if (dy > 0) r = 0.25;
327:       else r = 0.75;
328:     } else if (dx > 0) {
329:       r = PetscAtanReal((dy/ctx->vscale)/dx);
330:       if (dy >= 0) r /= 2*PETSC_PI;
331:       else r = r/(2*PETSC_PI)+1;
332:     } else r = PetscAtanReal((dy/ctx->vscale)/dx)/(2*PETSC_PI)+0.5;
333:     if (ctx->start_ang>ctx->end_ang) {
334:       if (r>ctx->end_ang && r<ctx->start_ang) *inside = -1;
335:     } else {
336:       if (r<ctx->start_ang || r>ctx->end_ang) *inside = -1;
337:     }
338:   }
339:   return(0);
340: }

342: PetscErrorCode RGSetFromOptions_Ring(PetscOptionItems *PetscOptionsObject,RG rg)
343: {
345:   PetscScalar    s;
346:   PetscReal      r1,r2,r3,r4,r5;
347:   PetscBool      flg1,flg2,flg3,flg4,flg5,flg6;

350:   PetscOptionsHead(PetscOptionsObject,"RG Ring Options");

352:     RGRingGetParameters(rg,&s,&r1,&r2,&r3,&r4,&r5);
353:     PetscOptionsScalar("-rg_ring_center","Center of ellipse","RGRingSetParameters",s,&s,&flg1);
354:     PetscOptionsReal("-rg_ring_radius","Radius of ellipse","RGRingSetParameters",r1,&r1,&flg2);
355:     PetscOptionsReal("-rg_ring_vscale","Vertical scale of ellipse","RGRingSetParameters",r2,&r2,&flg3);
356:     PetscOptionsReal("-rg_ring_startangle","Right-hand side angle","RGRingSetParameters",r3,&r3,&flg4);
357:     PetscOptionsReal("-rg_ring_endangle","Left-hand side angle","RGRingSetParameters",r4,&r4,&flg5);
358:     PetscOptionsReal("-rg_ring_width","Width of ring","RGRingSetParameters",r5,&r5,&flg6);
359:     if (flg1 || flg2 || flg3 || flg4 || flg5 || flg6) { RGRingSetParameters(rg,s,r1,r2,r3,r4,r5); }

361:   PetscOptionsTail();
362:   return(0);
363: }

365: PetscErrorCode RGDestroy_Ring(RG rg)
366: {

370:   PetscFree(rg->data);
371:   PetscObjectComposeFunction((PetscObject)rg,"RGRingSetParameters_C",NULL);
372:   PetscObjectComposeFunction((PetscObject)rg,"RGRingGetParameters_C",NULL);
373:   return(0);
374: }

376: SLEPC_EXTERN PetscErrorCode RGCreate_Ring(RG rg)
377: {
378:   RG_RING        *ring;

382:   PetscNewLog(rg,&ring);
383:   ring->center    = 0.0;
384:   ring->radius    = PETSC_MAX_REAL;
385:   ring->vscale    = 1.0;
386:   ring->start_ang = 0.0;
387:   ring->end_ang   = 1.0;
388:   ring->width     = 0.1;
389:   rg->data = (void*)ring;

391:   rg->ops->istrivial      = RGIsTrivial_Ring;
392:   rg->ops->computecontour = RGComputeContour_Ring;
393:   rg->ops->computebbox    = RGComputeBoundingBox_Ring;
394:   rg->ops->checkinside    = RGCheckInside_Ring;
395:   rg->ops->setfromoptions = RGSetFromOptions_Ring;
396:   rg->ops->view           = RGView_Ring;
397:   rg->ops->destroy        = RGDestroy_Ring;
398:   PetscObjectComposeFunction((PetscObject)rg,"RGRingSetParameters_C",RGRingSetParameters_Ring);
399:   PetscObjectComposeFunction((PetscObject)rg,"RGRingGetParameters_C",RGRingGetParameters_Ring);
400:   return(0);
401: }