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