Line data    Source code

       1             : /*
       2             :    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       3             :    SLEPc - Scalable Library for Eigenvalue Problem Computations
       4             :    Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
       5             : 
       6             :    This file is part of SLEPc.
       7             :    SLEPc is distributed under a 2-clause BSD license (see LICENSE).
       8             :    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       9             : */
      10             : /*
      11             :    Region enclosed in an ellipse (aligned with the coordinate axes)
      12             : */
      13             : 
      14             : #include <slepc/private/rgimpl.h>      /*I "slepcrg.h" I*/
      15             : #include <petscdraw.h>
      16             : 
      17             : typedef struct {
      18             :   PetscScalar center;     /* center of the ellipse */
      19             :   PetscReal   radius;     /* radius of the ellipse */
      20             :   PetscReal   vscale;     /* vertical scale of the ellipse */
      21             : } RG_ELLIPSE;
      22             : 
      23         140 : static PetscErrorCode RGEllipseSetParameters_Ellipse(RG rg,PetscScalar center,PetscReal radius,PetscReal vscale)
      24             : {
      25         140 :   RG_ELLIPSE *ctx = (RG_ELLIPSE*)rg->data;
      26             : 
      27         140 :   PetscFunctionBegin;
      28         140 :   ctx->center = center;
      29         140 :   if (radius == (PetscReal)PETSC_DETERMINE) {
      30           0 :     ctx->radius = 1.0;
      31         140 :   } else if (radius != (PetscReal)PETSC_CURRENT) {
      32         140 :     PetscCheck(radius>0.0,PetscObjectComm((PetscObject)rg),PETSC_ERR_ARG_OUTOFRANGE,"The radius argument must be > 0.0");
      33         140 :     ctx->radius = radius;
      34             :   }
      35         140 :   if (vscale == (PetscReal)PETSC_DETERMINE) {
      36           0 :     ctx->vscale = 1.0;
      37         140 :   } else if (vscale != (PetscReal)PETSC_CURRENT) {
      38         140 :     PetscCheck(vscale>0.0,PetscObjectComm((PetscObject)rg),PETSC_ERR_ARG_OUTOFRANGE,"The vscale argument must be > 0.0");
      39         140 :     ctx->vscale = vscale;
      40             :   }
      41         140 :   PetscFunctionReturn(PETSC_SUCCESS);
      42             : }
      43             : 
      44             : /*@
      45             :    RGEllipseSetParameters - Sets the parameters defining the ellipse region.
      46             : 
      47             :    Logically Collective
      48             : 
      49             :    Input Parameters:
      50             : +  rg     - the region context
      51             : .  center - center of the ellipse
      52             : .  radius - radius of the ellipse
      53             : -  vscale - vertical scale of the ellipse
      54             : 
      55             :    Options Database Keys:
      56             : +  -rg_ellipse_center - Sets the center
      57             : .  -rg_ellipse_radius - Sets the radius
      58             : -  -rg_ellipse_vscale - Sets the vertical scale
      59             : 
      60             :    Notes:
      61             :    In the case of complex scalars, a complex center can be provided in the
      62             :    command line with [+/-][realnumber][+/-]realnumberi with no spaces, e.g.
      63             :    -rg_ellipse_center 1.0+2.0i
      64             : 
      65             :    When PETSc is built with real scalars, the center is restricted to a real value.
      66             : 
      67             :    For radius and vscale, you can use PETSC_CURRENT to keep the current value, and
      68             :    PETSC_DETERMINE to set them to a default of 1.
      69             : 
      70             :    Level: advanced
      71             : 
      72             : .seealso: RGEllipseGetParameters()
      73             : @*/
      74         140 : PetscErrorCode RGEllipseSetParameters(RG rg,PetscScalar center,PetscReal radius,PetscReal vscale)
      75             : {
      76         140 :   PetscFunctionBegin;
      77         140 :   PetscValidHeaderSpecific(rg,RG_CLASSID,1);
      78         420 :   PetscValidLogicalCollectiveScalar(rg,center,2);
      79         420 :   PetscValidLogicalCollectiveReal(rg,radius,3);
      80         420 :   PetscValidLogicalCollectiveReal(rg,vscale,4);
      81         140 :   PetscTryMethod(rg,"RGEllipseSetParameters_C",(RG,PetscScalar,PetscReal,PetscReal),(rg,center,radius,vscale));
      82         140 :   PetscFunctionReturn(PETSC_SUCCESS);
      83             : }
      84             : 
      85         238 : static PetscErrorCode RGEllipseGetParameters_Ellipse(RG rg,PetscScalar *center,PetscReal *radius,PetscReal *vscale)
      86             : {
      87         238 :   RG_ELLIPSE *ctx = (RG_ELLIPSE*)rg->data;
      88             : 
      89         238 :   PetscFunctionBegin;
      90         238 :   if (center) *center = ctx->center;
      91         238 :   if (radius) *radius = ctx->radius;
      92         238 :   if (vscale) *vscale = ctx->vscale;
      93         238 :   PetscFunctionReturn(PETSC_SUCCESS);
      94             : }
      95             : 
      96             : /*@
      97             :    RGEllipseGetParameters - Gets the parameters that define the ellipse region.
      98             : 
      99             :    Not Collective
     100             : 
     101             :    Input Parameter:
     102             : .  rg     - the region context
     103             : 
     104             :    Output Parameters:
     105             : +  center - center of the region
     106             : .  radius - radius of the region
     107             : -  vscale - vertical scale of the region
     108             : 
     109             :    Level: advanced
     110             : 
     111             : .seealso: RGEllipseSetParameters()
     112             : @*/
     113         238 : PetscErrorCode RGEllipseGetParameters(RG rg,PetscScalar *center,PetscReal *radius,PetscReal *vscale)
     114             : {
     115         238 :   PetscFunctionBegin;
     116         238 :   PetscValidHeaderSpecific(rg,RG_CLASSID,1);
     117         238 :   PetscUseMethod(rg,"RGEllipseGetParameters_C",(RG,PetscScalar*,PetscReal*,PetscReal*),(rg,center,radius,vscale));
     118         238 :   PetscFunctionReturn(PETSC_SUCCESS);
     119             : }
     120             : 
     121           6 : static PetscErrorCode RGView_Ellipse(RG rg,PetscViewer viewer)
     122             : {
     123           6 :   RG_ELLIPSE     *ctx = (RG_ELLIPSE*)rg->data;
     124           6 :   PetscBool      isdraw,isascii;
     125           6 :   int            winw,winh;
     126           6 :   PetscDraw      draw;
     127           6 :   PetscDrawAxis  axis;
     128           6 :   PetscReal      cx,cy,r,ab,cd,lx,ly,w,scale=1.2;
     129           6 :   char           str[50];
     130             : 
     131           6 :   PetscFunctionBegin;
     132           6 :   PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw));
     133           6 :   PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
     134           6 :   if (isascii) {
     135           5 :     PetscCall(SlepcSNPrintfScalar(str,sizeof(str),ctx->center,PETSC_FALSE));
     136           5 :     PetscCall(PetscViewerASCIIPrintf(viewer,"  center: %s, radius: %g, vscale: %g\n",str,RGShowReal(ctx->radius),RGShowReal(ctx->vscale)));
     137           1 :   } else if (isdraw) {
     138           1 :     PetscCall(PetscViewerDrawGetDraw(viewer,0,&draw));
     139           1 :     PetscCall(PetscDrawCheckResizedWindow(draw));
     140           1 :     PetscCall(PetscDrawGetWindowSize(draw,&winw,&winh));
     141           1 :     winw = PetscMax(winw,1); winh = PetscMax(winh,1);
     142           1 :     PetscCall(PetscDrawClear(draw));
     143           1 :     PetscCall(PetscDrawSetTitle(draw,"Ellipse region"));
     144           1 :     PetscCall(PetscDrawAxisCreate(draw,&axis));
     145           1 :     cx = PetscRealPart(ctx->center)*rg->sfactor;
     146           1 :     cy = PetscImaginaryPart(ctx->center)*rg->sfactor;
     147           1 :     r  = ctx->radius*rg->sfactor;
     148           1 :     lx = 2*r;
     149           1 :     ly = 2*r*ctx->vscale;
     150           1 :     ab = cx;
     151           1 :     cd = cy;
     152           1 :     w  = scale*PetscMax(lx/winw,ly/winh)/2;
     153           1 :     PetscCall(PetscDrawAxisSetLimits(axis,ab-w*winw,ab+w*winw,cd-w*winh,cd+w*winh));
     154           1 :     PetscCall(PetscDrawAxisDraw(axis));
     155           1 :     PetscCall(PetscDrawAxisDestroy(&axis));
     156           1 :     PetscCall(PetscDrawEllipse(draw,cx,cy,2*r,2*r*ctx->vscale,PETSC_DRAW_RED));
     157           1 :     PetscCall(PetscDrawFlush(draw));
     158           1 :     PetscCall(PetscDrawSave(draw));
     159           1 :     PetscCall(PetscDrawPause(draw));
     160             :   }
     161           6 :   PetscFunctionReturn(PETSC_SUCCESS);
     162             : }
     163             : 
     164          83 : static PetscErrorCode RGIsTrivial_Ellipse(RG rg,PetscBool *trivial)
     165             : {
     166          83 :   RG_ELLIPSE *ctx = (RG_ELLIPSE*)rg->data;
     167             : 
     168          83 :   PetscFunctionBegin;
     169          83 :   if (rg->complement) *trivial = PetscNot(ctx->radius);
     170          83 :   else *trivial = PetscNot(ctx->radius<PETSC_MAX_REAL);
     171          83 :   PetscFunctionReturn(PETSC_SUCCESS);
     172             : }
     173             : 
     174         102 : static PetscErrorCode RGComputeContour_Ellipse(RG rg,PetscInt n,PetscScalar *cr,PetscScalar *ci)
     175             : {
     176         102 :   RG_ELLIPSE *ctx = (RG_ELLIPSE*)rg->data;
     177         102 :   PetscReal  theta;
     178         102 :   PetscInt   i;
     179             : 
     180         102 :   PetscFunctionBegin;
     181        4134 :   for (i=0;i<n;i++) {
     182        4032 :     theta = 2.0*PETSC_PI*(i+0.5)/n;
     183             : #if defined(PETSC_USE_COMPLEX)
     184        4032 :     cr[i] = ctx->center + ctx->radius*PetscCMPLX(PetscCosReal(theta),ctx->vscale*PetscSinReal(theta));
     185             : #else
     186             :     if (cr) cr[i] = ctx->center + ctx->radius*PetscCosReal(theta);
     187             :     if (ci) ci[i] = ctx->radius*ctx->vscale*PetscSinReal(theta);
     188             : #endif
     189             :   }
     190         102 :   PetscFunctionReturn(PETSC_SUCCESS);
     191             : }
     192             : 
     193          27 : static PetscErrorCode RGComputeBoundingBox_Ellipse(RG rg,PetscReal *a,PetscReal *b,PetscReal *c,PetscReal *d)
     194             : {
     195          27 :   RG_ELLIPSE *ctx = (RG_ELLIPSE*)rg->data;
     196             : 
     197          27 :   PetscFunctionBegin;
     198          27 :   if (a) *a = PetscRealPart(ctx->center) - ctx->radius;
     199          27 :   if (b) *b = PetscRealPart(ctx->center) + ctx->radius;
     200          27 :   if (c) *c = PetscImaginaryPart(ctx->center) - ctx->radius*ctx->vscale;
     201          27 :   if (d) *d = PetscImaginaryPart(ctx->center) + ctx->radius*ctx->vscale;
     202          27 :   PetscFunctionReturn(PETSC_SUCCESS);
     203             : }
     204             : 
     205         100 : static PetscErrorCode RGComputeQuadrature_Ellipse(RG rg,RGQuadRule quad,PetscInt n,PetscScalar *z,PetscScalar *zn,PetscScalar *w)
     206             : {
     207         100 :   RG_ELLIPSE *ctx = (RG_ELLIPSE*)rg->data;
     208         100 :   PetscReal  theta;
     209         100 :   PetscInt   i;
     210             : 
     211         100 :   PetscFunctionBegin;
     212        4112 :   for (i=0;i<n;i++) {
     213             : #if defined(PETSC_USE_COMPLEX)
     214        4012 :     theta = 2.0*PETSC_PI*(i+0.5)/n;
     215        4012 :     zn[i] = PetscCMPLX(PetscCosReal(theta),ctx->vscale*PetscSinReal(theta));
     216        4012 :     w[i]  = rg->sfactor*ctx->radius*(PetscCMPLX(ctx->vscale*PetscCosReal(theta),PetscSinReal(theta)))/n;
     217             : #else
     218             :     theta = PETSC_PI*(i+0.5)/n;
     219             :     zn[i] = PetscCosReal(theta);
     220             :     w[i]  = PetscCosReal((n-1)*theta)/n;
     221             :     z[i]  = rg->sfactor*(ctx->center + ctx->radius*zn[i]);
     222             : #endif
     223             :   }
     224         100 :   PetscFunctionReturn(PETSC_SUCCESS);
     225             : }
     226             : 
     227        3640 : static PetscErrorCode RGCheckInside_Ellipse(RG rg,PetscReal px,PetscReal py,PetscInt *inside)
     228             : {
     229        3640 :   RG_ELLIPSE *ctx = (RG_ELLIPSE*)rg->data;
     230        3640 :   PetscReal  dx,dy,r;
     231             : 
     232        3640 :   PetscFunctionBegin;
     233             : #if defined(PETSC_USE_COMPLEX)
     234        3640 :   dx = (px-PetscRealPart(ctx->center))/ctx->radius;
     235        3640 :   dy = (py-PetscImaginaryPart(ctx->center))/ctx->radius;
     236             : #else
     237             :   dx = (px-ctx->center)/ctx->radius;
     238             :   dy = py/ctx->radius;
     239             : #endif
     240        3640 :   r = 1.0-dx*dx-(dy*dy)/(ctx->vscale*ctx->vscale);
     241        3640 :   *inside = PetscSign(r);
     242        3640 :   PetscFunctionReturn(PETSC_SUCCESS);
     243             : }
     244             : 
     245          30 : static PetscErrorCode RGIsAxisymmetric_Ellipse(RG rg,PetscBool vertical,PetscBool *symm)
     246             : {
     247          30 :   RG_ELLIPSE *ctx = (RG_ELLIPSE*)rg->data;
     248             : 
     249          30 :   PetscFunctionBegin;
     250          30 :   if (vertical) *symm = (PetscRealPart(ctx->center) == 0.0)? PETSC_TRUE: PETSC_FALSE;
     251          28 :   else *symm = (PetscImaginaryPart(ctx->center) == 0.0)? PETSC_TRUE: PETSC_FALSE;
     252          30 :   PetscFunctionReturn(PETSC_SUCCESS);
     253             : }
     254             : 
     255         152 : static PetscErrorCode RGSetFromOptions_Ellipse(RG rg,PetscOptionItems *PetscOptionsObject)
     256             : {
     257         152 :   PetscScalar    s;
     258         152 :   PetscReal      r1,r2;
     259         152 :   PetscBool      flg1,flg2,flg3;
     260             : 
     261         152 :   PetscFunctionBegin;
     262         152 :   PetscOptionsHeadBegin(PetscOptionsObject,"RG Ellipse Options");
     263             : 
     264         152 :     PetscCall(RGEllipseGetParameters(rg,&s,&r1,&r2));
     265         152 :     PetscCall(PetscOptionsScalar("-rg_ellipse_center","Center of ellipse","RGEllipseSetParameters",s,&s,&flg1));
     266         152 :     PetscCall(PetscOptionsReal("-rg_ellipse_radius","Radius of ellipse","RGEllipseSetParameters",r1,&r1,&flg2));
     267         152 :     PetscCall(PetscOptionsReal("-rg_ellipse_vscale","Vertical scale of ellipse","RGEllipseSetParameters",r2,&r2,&flg3));
     268         152 :     if (flg1 || flg2 || flg3) PetscCall(RGEllipseSetParameters(rg,s,r1,r2));
     269             : 
     270         152 :   PetscOptionsHeadEnd();
     271         152 :   PetscFunctionReturn(PETSC_SUCCESS);
     272             : }
     273             : 
     274          83 : static PetscErrorCode RGDestroy_Ellipse(RG rg)
     275             : {
     276          83 :   PetscFunctionBegin;
     277          83 :   PetscCall(PetscFree(rg->data));
     278          83 :   PetscCall(PetscObjectComposeFunction((PetscObject)rg,"RGEllipseSetParameters_C",NULL));
     279          83 :   PetscCall(PetscObjectComposeFunction((PetscObject)rg,"RGEllipseGetParameters_C",NULL));
     280          83 :   PetscFunctionReturn(PETSC_SUCCESS);
     281             : }
     282             : 
     283          83 : SLEPC_EXTERN PetscErrorCode RGCreate_Ellipse(RG rg)
     284             : {
     285          83 :   RG_ELLIPSE     *ellipse;
     286             : 
     287          83 :   PetscFunctionBegin;
     288          83 :   PetscCall(PetscNew(&ellipse));
     289          83 :   ellipse->center = 0.0;
     290          83 :   ellipse->radius = PETSC_MAX_REAL;
     291          83 :   ellipse->vscale = 1.0;
     292          83 :   rg->data = (void*)ellipse;
     293             : 
     294          83 :   rg->ops->istrivial         = RGIsTrivial_Ellipse;
     295          83 :   rg->ops->computecontour    = RGComputeContour_Ellipse;
     296          83 :   rg->ops->computebbox       = RGComputeBoundingBox_Ellipse;
     297          83 :   rg->ops->computequadrature = RGComputeQuadrature_Ellipse;
     298          83 :   rg->ops->checkinside       = RGCheckInside_Ellipse;
     299          83 :   rg->ops->isaxisymmetric    = RGIsAxisymmetric_Ellipse;
     300          83 :   rg->ops->setfromoptions    = RGSetFromOptions_Ellipse;
     301          83 :   rg->ops->view              = RGView_Ellipse;
     302          83 :   rg->ops->destroy           = RGDestroy_Ellipse;
     303          83 :   PetscCall(PetscObjectComposeFunction((PetscObject)rg,"RGEllipseSetParameters_C",RGEllipseSetParameters_Ellipse));
     304          83 :   PetscCall(PetscObjectComposeFunction((PetscObject)rg,"RGEllipseGetParameters_C",RGEllipseGetParameters_Ellipse));
     305          83 :   PetscFunctionReturn(PETSC_SUCCESS);
     306             : }