Actual source code: rgring.c

slepc-3.18.1 2022-11-02
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-, 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>
 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;

 31:   ctx->center = center;
 32:   if (radius == PETSC_DEFAULT) {
 33:     ctx->radius = 1.0;
 34:   } else {
 36:     ctx->radius = radius;
 37:   }
 39:   ctx->vscale = vscale;
 40:   if (start_ang == PETSC_DEFAULT) {
 41:     ctx->start_ang = 0.0;
 42:   } else {
 45:     ctx->start_ang = start_ang;
 46:   }
 47:   if (end_ang == PETSC_DEFAULT) {
 48:     ctx->end_ang = 1.0;
 49:   } else {
 52:     ctx->end_ang = end_ang;
 53:   }
 54: #if !defined(PETSC_USE_COMPLEX)
 56: #endif
 57:   if (width == PETSC_DEFAULT) {
 58:     ctx->width = 0.1;
 59:   } else {
 61:     ctx->width = width;
 62:   }
 63:   return 0;
 64: }

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

 69:    Logically Collective on rg

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

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

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

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

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

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

108:    Level: advanced

110: .seealso: RGRingGetParameters()
111: @*/
112: PetscErrorCode RGRingSetParameters(RG rg,PetscScalar center,PetscReal radius,PetscReal vscale,PetscReal start_ang,PetscReal end_ang,PetscReal width)
113: {
121:   PetscTryMethod(rg,"RGRingSetParameters_C",(RG,PetscScalar,PetscReal,PetscReal,PetscReal,PetscReal,PetscReal),(rg,center,radius,vscale,start_ang,end_ang,width));
122:   return 0;
123: }

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

129:   if (center)    *center    = ctx->center;
130:   if (radius)    *radius    = ctx->radius;
131:   if (vscale)    *vscale    = ctx->vscale;
132:   if (start_ang) *start_ang = ctx->start_ang;
133:   if (end_ang)   *end_ang   = ctx->end_ang;
134:   if (width)     *width     = ctx->width;
135:   return 0;
136: }

138: /*@
139:    RGRingGetParameters - Gets the parameters that define the ring region.

141:    Not Collective

143:    Input Parameter:
144: .  rg     - the region context

146:    Output Parameters:
147: +  center    - center of the region
148: .  radius    - radius of the region
149: .  vscale    - vertical scale of the region
150: .  start_ang - the right-hand side angle
151: .  end_ang   - the left-hand side angle
152: -  width     - width of the ring

154:    Level: advanced

156: .seealso: RGRingSetParameters()
157: @*/
158: PetscErrorCode RGRingGetParameters(RG rg,PetscScalar *center,PetscReal *radius,PetscReal *vscale,PetscReal *start_ang,PetscReal *end_ang,PetscReal *width)
159: {
161:   PetscUseMethod(rg,"RGRingGetParameters_C",(RG,PetscScalar*,PetscReal*,PetscReal*,PetscReal*,PetscReal*,PetscReal*),(rg,center,radius,vscale,start_ang,end_ang,width));
162:   return 0;
163: }

165: PetscErrorCode RGView_Ring(RG rg,PetscViewer viewer)
166: {
167:   RG_RING        *ctx = (RG_RING*)rg->data;
168:   int            winw,winh;
169:   PetscBool      isdraw,isascii;
170:   PetscDraw      draw;
171:   PetscDrawAxis  axis;
172:   PetscReal      cx,cy,radius,width,ab,cd,lx,ly,w,end_ang,x1,y1,x2,y2,r,theta,scale=1.2;
173:   char           str[50];

175:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
176:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
177:   if (isascii) {
178:     SlepcSNPrintfScalar(str,sizeof(str),ctx->center,PETSC_FALSE);
179:     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);
180:   } else if (isdraw) {
181:     PetscViewerDrawGetDraw(viewer,0,&draw);
182:     PetscDrawCheckResizedWindow(draw);
183:     PetscDrawGetWindowSize(draw,&winw,&winh);
184:     winw = PetscMax(winw,1); winh = PetscMax(winh,1);
185:     PetscDrawClear(draw);
186:     PetscDrawSetTitle(draw,"Ring region");
187:     PetscDrawAxisCreate(draw,&axis);
188:     cx = PetscRealPart(ctx->center)*rg->sfactor;
189:     cy = PetscImaginaryPart(ctx->center)*rg->sfactor;
190:     radius = ctx->radius*rg->sfactor;
191:     width  = ctx->width*rg->sfactor;
192:     lx = 2*(radius+width);
193:     ly = 2*(radius+width)*ctx->vscale;
194:     ab = cx;
195:     cd = cy;
196:     w  = scale*PetscMax(lx/winw,ly/winh)/2;
197:     PetscDrawAxisSetLimits(axis,ab-w*winw,ab+w*winw,cd-w*winh,cd+w*winh);
198:     PetscDrawAxisDraw(axis);
199:     PetscDrawAxisDestroy(&axis);
200:     /* draw outer ellipse */
201:     PetscDrawEllipse(draw,cx,cy,2*(radius+width),2*(radius+width)*ctx->vscale,PETSC_DRAW_ORANGE);
202:     /* remove inner part */
203:     PetscDrawEllipse(draw,cx,cy,2*(radius-width),2*(radius-width)*ctx->vscale,PETSC_DRAW_WHITE);
204:     if (ctx->start_ang!=ctx->end_ang) {
205:       /* remove section from end_ang to start_ang */
206:       end_ang = (ctx->start_ang<ctx->end_ang)? ctx->end_ang-1: ctx->end_ang;
207:       theta = end_ang;
208:       r = scale*(radius+width);
209:       if (ctx->vscale>1) r *= ctx->vscale;
210:       x1 = PetscMin(PetscMax(ab+r*PetscCosReal(2.0*PETSC_PI*theta),ab-w*winw),ab+w*winw);
211:       y1 = PetscMin(PetscMax(cd+r*PetscSinReal(2.0*PETSC_PI*theta),cd-w*winh),cd+w*winh);
212:       do {
213:         theta = PetscMin(PetscFloorReal(8*theta+1)/8,ctx->start_ang);
214:         x2 = PetscMin(PetscMax(ab+r*PetscCosReal(2.0*PETSC_PI*theta),ab-w*winw),ab+w*winw);
215:         y2 = PetscMin(PetscMax(cd+r*PetscSinReal(2.0*PETSC_PI*theta),cd-w*winh),cd+w*winh);
216:         PetscDrawTriangle(draw,cx,cy,x1,y1,x2,y2,PETSC_DRAW_WHITE,PETSC_DRAW_WHITE,PETSC_DRAW_WHITE);
217:         x1 = x2; y1 = y2;
218:       } while (theta<ctx->start_ang);
219:     }
220:     PetscDrawFlush(draw);
221:     PetscDrawSave(draw);
222:     PetscDrawPause(draw);
223:   }
224:   return 0;
225: }

227: PetscErrorCode RGIsTrivial_Ring(RG rg,PetscBool *trivial)
228: {
229:   RG_RING *ctx = (RG_RING*)rg->data;

231:   if (rg->complement) *trivial = PetscNot(ctx->radius);
232:   else *trivial = PetscNot(ctx->radius<PETSC_MAX_REAL);
233:   return 0;
234: }

236: PetscErrorCode RGComputeContour_Ring(RG rg,PetscInt n,PetscScalar *cr,PetscScalar *ci)
237: {
238:   RG_RING   *ctx = (RG_RING*)rg->data;
239:   PetscReal theta,start_ang;
240:   PetscInt  i,n2=n/2;

242:   start_ang = (ctx->start_ang>ctx->end_ang)? ctx->start_ang-1: ctx->start_ang;
243:   for (i=0;i<n;i++) {
244:     if (i < n2) {
245:       theta = ((ctx->end_ang-start_ang)*i/n2 + start_ang)*2.0*PETSC_PI;
246: #if defined(PETSC_USE_COMPLEX)
247:       cr[i] = ctx->center + (ctx->radius+ctx->width/2.0)*PetscCMPLX(PetscCosReal(theta),ctx->vscale*PetscSinReal(theta));
248: #else
249:       if (cr) cr[i] = ctx->center + (ctx->radius+ctx->width/2.0)*PetscCosReal(theta);
250:       if (ci) ci[i] = (ctx->radius+ctx->width/2.0)*ctx->vscale*PetscSinReal(theta);
251: #endif
252:     } else {
253:       theta = ((ctx->end_ang-start_ang)*(n-i)/n2 + start_ang)*2.0*PETSC_PI;
254: #if defined(PETSC_USE_COMPLEX)
255:       cr[i] = ctx->center + (ctx->radius-ctx->width/2.0)*PetscCMPLX(PetscCosReal(theta),ctx->vscale*PetscSinReal(theta));
256: #else
257:       if (cr) cr[i] = ctx->center + (ctx->radius-ctx->width/2.0)*PetscCosReal(theta);
258:       if (ci) ci[i] = (ctx->radius-ctx->width/2.0)*ctx->vscale*PetscSinReal(theta);
259: #endif
260:     }
261:   }
262:   return 0;
263: }

265: PetscErrorCode RGComputeBoundingBox_Ring(RG rg,PetscReal *a,PetscReal *b,PetscReal *c,PetscReal *d)
266: {
267:   RG_RING *ctx = (RG_RING*)rg->data;

269:   /* current implementation does not return a tight bounding box */
270:   if (a) *a = PetscRealPart(ctx->center) - (ctx->radius+ctx->width/2.0);
271:   if (b) *b = PetscRealPart(ctx->center) + (ctx->radius+ctx->width/2.0);
272:   if (c) *c = PetscImaginaryPart(ctx->center) - (ctx->radius+ctx->width/2.0)*ctx->vscale;
273:   if (d) *d = PetscImaginaryPart(ctx->center) + (ctx->radius+ctx->width/2.0)*ctx->vscale;
274:   return 0;
275: }

277: PetscErrorCode RGComputeQuadrature_Ring(RG rg,RGQuadRule quad,PetscInt n,PetscScalar *z,PetscScalar *zn,PetscScalar *w)
278: {
279:   RG_RING     *ctx = (RG_RING*)rg->data;
280:   PetscReal   max_w=0.0;
281:   PetscScalar tmp,tmp2;
282:   PetscInt    i,j;

284:   if (quad == RG_QUADRULE_CHEBYSHEV) {
285: #if defined(PETSC_USE_COMPLEX)
286:     PetscReal theta;
287:     for (i=0;i<n;i++) {
288:       theta = PETSC_PI*(i+0.5)/n;
289:       zn[i] = PetscCosReal(theta);
290:       w[i]  = PetscCosReal((n-1)*theta)/n;
291:       theta = (ctx->start_ang*2.0+(ctx->end_ang-ctx->start_ang)*(PetscRealPart(zn[i])+1.0))*PETSC_PI;
292:       z[i] = rg->sfactor*(ctx->center + ctx->radius*PetscCMPLX(PetscCosReal(theta),ctx->vscale*PetscSinReal(theta)));
293:     }
294: #else
295:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Integration points on a vertical line require complex arithmetic");
296: #endif
297:   } else {  /* RG_QUADRULE_TRAPEZOIDAL */
298:     for (i=0;i<n;i++) {
299:       zn[i] = (z[i]-rg->sfactor*ctx->center)/(rg->sfactor*ctx->radius);
300:       tmp = 1.0; tmp2 = 1.0;
301:       for (j=0;j<n;j++) {
302:         tmp *= z[j];
303:         if (i != j) tmp2 *= z[j]-z[i];
304:       }
305:       w[i] = tmp/tmp2;
306:       max_w = PetscMax(PetscAbsScalar(w[i]),max_w);
307:     }
308:     for (i=0;i<n;i++) w[i] /= (PetscScalar)max_w;
309:   }
310:   return 0;
311: }

313: PetscErrorCode RGCheckInside_Ring(RG rg,PetscReal px,PetscReal py,PetscInt *inside)
314: {
315:   RG_RING   *ctx = (RG_RING*)rg->data;
316:   PetscReal dx,dy,r;

318:   /* outer ellipse */
319: #if defined(PETSC_USE_COMPLEX)
320:   dx = (px-PetscRealPart(ctx->center))/(ctx->radius+ctx->width/2.0);
321:   dy = (py-PetscImaginaryPart(ctx->center))/(ctx->radius+ctx->width/2.0);
322: #else
323:   dx = (px-ctx->center)/(ctx->radius+ctx->width/2.0);
324:   dy = py/(ctx->radius+ctx->width/2.0);
325: #endif
326:   r = 1.0-dx*dx-(dy*dy)/(ctx->vscale*ctx->vscale);
327:   *inside = PetscSign(r);
328:   /* inner ellipse */
329: #if defined(PETSC_USE_COMPLEX)
330:   dx = (px-PetscRealPart(ctx->center))/(ctx->radius-ctx->width/2.0);
331:   dy = (py-PetscImaginaryPart(ctx->center))/(ctx->radius-ctx->width/2.0);
332: #else
333:   dx = (px-ctx->center)/(ctx->radius-ctx->width/2.0);
334:   dy = py/(ctx->radius-ctx->width/2.0);
335: #endif
336:   r = -1.0+dx*dx+(dy*dy)/(ctx->vscale*ctx->vscale);
337:   *inside *= PetscSign(r);
338:   if (*inside == 1) {  /* check angles */
339: #if defined(PETSC_USE_COMPLEX)
340:     dx = (px-PetscRealPart(ctx->center));
341:     dy = (py-PetscImaginaryPart(ctx->center));
342: #else
343:     dx = px-ctx->center;
344:     dy = py;
345: #endif
346:     if (dx == 0) {
347:       if (dy == 0) r = -1;
348:       else if (dy > 0) r = 0.25;
349:       else r = 0.75;
350:     } else if (dx > 0) {
351:       r = PetscAtanReal((dy/ctx->vscale)/dx);
352:       if (dy >= 0) r /= 2*PETSC_PI;
353:       else r = r/(2*PETSC_PI)+1;
354:     } else r = PetscAtanReal((dy/ctx->vscale)/dx)/(2*PETSC_PI)+0.5;
355:     if (ctx->start_ang>ctx->end_ang) {
356:       if (r>ctx->end_ang && r<ctx->start_ang) *inside = -1;
357:     } else {
358:       if (r<ctx->start_ang || r>ctx->end_ang) *inside = -1;
359:     }
360:   }
361:   return 0;
362: }

364: PetscErrorCode RGIsAxisymmetric_Ring(RG rg,PetscBool vertical,PetscBool *symm)
365: {
366:   RG_RING *ctx = (RG_RING*)rg->data;

368:   if (vertical) *symm = (PetscRealPart(ctx->center) == 0.0 && PetscAbs(ctx->start_ang+ctx->end_ang-PetscRealConstant(1.0)) == 0.5)? PETSC_TRUE: PETSC_FALSE;
369:   else *symm = (PetscImaginaryPart(ctx->center) == 0.0 && ctx->start_ang+ctx->end_ang == 1.0)? PETSC_TRUE: PETSC_FALSE;
370:   return 0;
371: }

373: PetscErrorCode RGSetFromOptions_Ring(RG rg,PetscOptionItems *PetscOptionsObject)
374: {
375:   PetscScalar    s;
376:   PetscReal      r1,r2,r3,r4,r5;
377:   PetscBool      flg1,flg2,flg3,flg4,flg5,flg6;

379:   PetscOptionsHeadBegin(PetscOptionsObject,"RG Ring Options");

381:     RGRingGetParameters(rg,&s,&r1,&r2,&r3,&r4,&r5);
382:     PetscOptionsScalar("-rg_ring_center","Center of ellipse","RGRingSetParameters",s,&s,&flg1);
383:     PetscOptionsReal("-rg_ring_radius","Radius of ellipse","RGRingSetParameters",r1,&r1,&flg2);
384:     PetscOptionsReal("-rg_ring_vscale","Vertical scale of ellipse","RGRingSetParameters",r2,&r2,&flg3);
385:     PetscOptionsReal("-rg_ring_startangle","Right-hand side angle","RGRingSetParameters",r3,&r3,&flg4);
386:     PetscOptionsReal("-rg_ring_endangle","Left-hand side angle","RGRingSetParameters",r4,&r4,&flg5);
387:     PetscOptionsReal("-rg_ring_width","Width of ring","RGRingSetParameters",r5,&r5,&flg6);
388:     if (flg1 || flg2 || flg3 || flg4 || flg5 || flg6) RGRingSetParameters(rg,s,r1,r2,r3,r4,r5);

390:   PetscOptionsHeadEnd();
391:   return 0;
392: }

394: PetscErrorCode RGDestroy_Ring(RG rg)
395: {
396:   PetscFree(rg->data);
397:   PetscObjectComposeFunction((PetscObject)rg,"RGRingSetParameters_C",NULL);
398:   PetscObjectComposeFunction((PetscObject)rg,"RGRingGetParameters_C",NULL);
399:   return 0;
400: }

402: SLEPC_EXTERN PetscErrorCode RGCreate_Ring(RG rg)
403: {
404:   RG_RING        *ring;

406:   PetscNew(&ring);
407:   ring->center    = 0.0;
408:   ring->radius    = PETSC_MAX_REAL;
409:   ring->vscale    = 1.0;
410:   ring->start_ang = 0.0;
411:   ring->end_ang   = 1.0;
412:   ring->width     = 0.1;
413:   rg->data = (void*)ring;

415:   rg->ops->istrivial         = RGIsTrivial_Ring;
416:   rg->ops->computecontour    = RGComputeContour_Ring;
417:   rg->ops->computebbox       = RGComputeBoundingBox_Ring;
418:   rg->ops->computequadrature = RGComputeQuadrature_Ring;
419:   rg->ops->checkinside       = RGCheckInside_Ring;
420:   rg->ops->isaxisymmetric    = RGIsAxisymmetric_Ring;
421:   rg->ops->setfromoptions    = RGSetFromOptions_Ring;
422:   rg->ops->view              = RGView_Ring;
423:   rg->ops->destroy           = RGDestroy_Ring;
424:   PetscObjectComposeFunction((PetscObject)rg,"RGRingSetParameters_C",RGRingSetParameters_Ring);
425:   PetscObjectComposeFunction((PetscObject)rg,"RGRingGetParameters_C",RGRingGetParameters_Ring);
426:   return 0;
427: }