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 41 : static PetscErrorCode RGEllipseSetParameters_Ellipse(RG rg,PetscScalar center,PetscReal radius,PetscReal vscale)
24 : {
25 41 : RG_ELLIPSE *ctx = (RG_ELLIPSE*)rg->data;
26 :
27 41 : PetscFunctionBegin;
28 41 : ctx->center = center;
29 41 : if (radius == (PetscReal)PETSC_DETERMINE) {
30 0 : ctx->radius = 1.0;
31 41 : } else if (radius != (PetscReal)PETSC_CURRENT) {
32 41 : PetscCheck(radius>0.0,PetscObjectComm((PetscObject)rg),PETSC_ERR_ARG_OUTOFRANGE,"The radius argument must be > 0.0");
33 41 : ctx->radius = radius;
34 : }
35 41 : if (vscale == (PetscReal)PETSC_DETERMINE) {
36 0 : ctx->vscale = 1.0;
37 41 : } else if (vscale != (PetscReal)PETSC_CURRENT) {
38 41 : PetscCheck(vscale>0.0,PetscObjectComm((PetscObject)rg),PETSC_ERR_ARG_OUTOFRANGE,"The vscale argument must be > 0.0");
39 41 : ctx->vscale = vscale;
40 : }
41 41 : 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 41 : PetscErrorCode RGEllipseSetParameters(RG rg,PetscScalar center,PetscReal radius,PetscReal vscale)
75 : {
76 41 : PetscFunctionBegin;
77 41 : PetscValidHeaderSpecific(rg,RG_CLASSID,1);
78 123 : PetscValidLogicalCollectiveScalar(rg,center,2);
79 123 : PetscValidLogicalCollectiveReal(rg,radius,3);
80 123 : PetscValidLogicalCollectiveReal(rg,vscale,4);
81 41 : PetscTryMethod(rg,"RGEllipseSetParameters_C",(RG,PetscScalar,PetscReal,PetscReal),(rg,center,radius,vscale));
82 41 : PetscFunctionReturn(PETSC_SUCCESS);
83 : }
84 :
85 44 : static PetscErrorCode RGEllipseGetParameters_Ellipse(RG rg,PetscScalar *center,PetscReal *radius,PetscReal *vscale)
86 : {
87 44 : RG_ELLIPSE *ctx = (RG_ELLIPSE*)rg->data;
88 :
89 44 : PetscFunctionBegin;
90 44 : if (center) *center = ctx->center;
91 44 : if (radius) *radius = ctx->radius;
92 44 : if (vscale) *vscale = ctx->vscale;
93 44 : 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 44 : PetscErrorCode RGEllipseGetParameters(RG rg,PetscScalar *center,PetscReal *radius,PetscReal *vscale)
114 : {
115 44 : PetscFunctionBegin;
116 44 : PetscValidHeaderSpecific(rg,RG_CLASSID,1);
117 44 : PetscUseMethod(rg,"RGEllipseGetParameters_C",(RG,PetscScalar*,PetscReal*,PetscReal*),(rg,center,radius,vscale));
118 44 : PetscFunctionReturn(PETSC_SUCCESS);
119 : }
120 :
121 3 : static PetscErrorCode RGView_Ellipse(RG rg,PetscViewer viewer)
122 : {
123 3 : RG_ELLIPSE *ctx = (RG_ELLIPSE*)rg->data;
124 3 : PetscBool isdraw,isascii;
125 3 : int winw,winh;
126 3 : PetscDraw draw;
127 3 : PetscDrawAxis axis;
128 3 : PetscReal cx,cy,r,ab,cd,lx,ly,w,scale=1.2;
129 3 : char str[50];
130 :
131 3 : PetscFunctionBegin;
132 3 : PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw));
133 3 : PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
134 3 : if (isascii) {
135 2 : PetscCall(SlepcSNPrintfScalar(str,sizeof(str),ctx->center,PETSC_FALSE));
136 2 : 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 3 : PetscFunctionReturn(PETSC_SUCCESS);
162 : }
163 :
164 25 : static PetscErrorCode RGIsTrivial_Ellipse(RG rg,PetscBool *trivial)
165 : {
166 25 : RG_ELLIPSE *ctx = (RG_ELLIPSE*)rg->data;
167 :
168 25 : PetscFunctionBegin;
169 25 : if (rg->complement) *trivial = PetscNot(ctx->radius);
170 25 : else *trivial = PetscNot(ctx->radius<PETSC_MAX_REAL);
171 25 : PetscFunctionReturn(PETSC_SUCCESS);
172 : }
173 :
174 19 : static PetscErrorCode RGComputeContour_Ellipse(RG rg,PetscInt n,PetscScalar *cr,PetscScalar *ci)
175 : {
176 19 : RG_ELLIPSE *ctx = (RG_ELLIPSE*)rg->data;
177 19 : PetscReal theta;
178 19 : PetscInt i;
179 :
180 19 : PetscFunctionBegin;
181 559 : for (i=0;i<n;i++) {
182 540 : theta = 2.0*PETSC_PI*(i+0.5)/n;
183 : #if defined(PETSC_USE_COMPLEX)
184 : cr[i] = ctx->center + ctx->radius*PetscCMPLX(PetscCosReal(theta),ctx->vscale*PetscSinReal(theta));
185 : #else
186 540 : if (cr) cr[i] = ctx->center + ctx->radius*PetscCosReal(theta);
187 540 : if (ci) ci[i] = ctx->radius*ctx->vscale*PetscSinReal(theta);
188 : #endif
189 : }
190 19 : PetscFunctionReturn(PETSC_SUCCESS);
191 : }
192 :
193 2 : static PetscErrorCode RGComputeBoundingBox_Ellipse(RG rg,PetscReal *a,PetscReal *b,PetscReal *c,PetscReal *d)
194 : {
195 2 : RG_ELLIPSE *ctx = (RG_ELLIPSE*)rg->data;
196 :
197 2 : PetscFunctionBegin;
198 2 : if (a) *a = PetscRealPart(ctx->center) - ctx->radius;
199 2 : if (b) *b = PetscRealPart(ctx->center) + ctx->radius;
200 2 : if (c) *c = PetscImaginaryPart(ctx->center) - ctx->radius*ctx->vscale;
201 2 : if (d) *d = PetscImaginaryPart(ctx->center) + ctx->radius*ctx->vscale;
202 2 : PetscFunctionReturn(PETSC_SUCCESS);
203 : }
204 :
205 17 : static PetscErrorCode RGComputeQuadrature_Ellipse(RG rg,RGQuadRule quad,PetscInt n,PetscScalar *z,PetscScalar *zn,PetscScalar *w)
206 : {
207 17 : RG_ELLIPSE *ctx = (RG_ELLIPSE*)rg->data;
208 17 : PetscReal theta;
209 17 : PetscInt i;
210 :
211 17 : PetscFunctionBegin;
212 537 : for (i=0;i<n;i++) {
213 : #if defined(PETSC_USE_COMPLEX)
214 : theta = 2.0*PETSC_PI*(i+0.5)/n;
215 : zn[i] = PetscCMPLX(PetscCosReal(theta),ctx->vscale*PetscSinReal(theta));
216 : w[i] = rg->sfactor*ctx->radius*(PetscCMPLX(ctx->vscale*PetscCosReal(theta),PetscSinReal(theta)))/n;
217 : #else
218 520 : theta = PETSC_PI*(i+0.5)/n;
219 520 : zn[i] = PetscCosReal(theta);
220 520 : w[i] = PetscCosReal((n-1)*theta)/n;
221 520 : z[i] = rg->sfactor*(ctx->center + ctx->radius*zn[i]);
222 : #endif
223 : }
224 17 : PetscFunctionReturn(PETSC_SUCCESS);
225 : }
226 :
227 1342 : static PetscErrorCode RGCheckInside_Ellipse(RG rg,PetscReal px,PetscReal py,PetscInt *inside)
228 : {
229 1342 : RG_ELLIPSE *ctx = (RG_ELLIPSE*)rg->data;
230 1342 : PetscReal dx,dy,r;
231 :
232 1342 : PetscFunctionBegin;
233 : #if defined(PETSC_USE_COMPLEX)
234 : dx = (px-PetscRealPart(ctx->center))/ctx->radius;
235 : dy = (py-PetscImaginaryPart(ctx->center))/ctx->radius;
236 : #else
237 1342 : dx = (px-ctx->center)/ctx->radius;
238 1342 : dy = py/ctx->radius;
239 : #endif
240 1342 : r = 1.0-dx*dx-(dy*dy)/(ctx->vscale*ctx->vscale);
241 1342 : *inside = PetscSign(r);
242 1342 : PetscFunctionReturn(PETSC_SUCCESS);
243 : }
244 :
245 2 : static PetscErrorCode RGIsAxisymmetric_Ellipse(RG rg,PetscBool vertical,PetscBool *symm)
246 : {
247 2 : RG_ELLIPSE *ctx = (RG_ELLIPSE*)rg->data;
248 :
249 2 : PetscFunctionBegin;
250 2 : if (vertical) *symm = (PetscRealPart(ctx->center) == 0.0)? PETSC_TRUE: PETSC_FALSE;
251 1 : else *symm = (PetscImaginaryPart(ctx->center) == 0.0)? PETSC_TRUE: PETSC_FALSE;
252 2 : PetscFunctionReturn(PETSC_SUCCESS);
253 : }
254 :
255 40 : static PetscErrorCode RGSetFromOptions_Ellipse(RG rg,PetscOptionItems *PetscOptionsObject)
256 : {
257 40 : PetscScalar s;
258 40 : PetscReal r1,r2;
259 40 : PetscBool flg1,flg2,flg3;
260 :
261 40 : PetscFunctionBegin;
262 40 : PetscOptionsHeadBegin(PetscOptionsObject,"RG Ellipse Options");
263 :
264 40 : PetscCall(RGEllipseGetParameters(rg,&s,&r1,&r2));
265 40 : PetscCall(PetscOptionsScalar("-rg_ellipse_center","Center of ellipse","RGEllipseSetParameters",s,&s,&flg1));
266 40 : PetscCall(PetscOptionsReal("-rg_ellipse_radius","Radius of ellipse","RGEllipseSetParameters",r1,&r1,&flg2));
267 40 : PetscCall(PetscOptionsReal("-rg_ellipse_vscale","Vertical scale of ellipse","RGEllipseSetParameters",r2,&r2,&flg3));
268 40 : if (flg1 || flg2 || flg3) PetscCall(RGEllipseSetParameters(rg,s,r1,r2));
269 :
270 40 : PetscOptionsHeadEnd();
271 40 : PetscFunctionReturn(PETSC_SUCCESS);
272 : }
273 :
274 25 : static PetscErrorCode RGDestroy_Ellipse(RG rg)
275 : {
276 25 : PetscFunctionBegin;
277 25 : PetscCall(PetscFree(rg->data));
278 25 : PetscCall(PetscObjectComposeFunction((PetscObject)rg,"RGEllipseSetParameters_C",NULL));
279 25 : PetscCall(PetscObjectComposeFunction((PetscObject)rg,"RGEllipseGetParameters_C",NULL));
280 25 : PetscFunctionReturn(PETSC_SUCCESS);
281 : }
282 :
283 25 : SLEPC_EXTERN PetscErrorCode RGCreate_Ellipse(RG rg)
284 : {
285 25 : RG_ELLIPSE *ellipse;
286 :
287 25 : PetscFunctionBegin;
288 25 : PetscCall(PetscNew(&ellipse));
289 25 : ellipse->center = 0.0;
290 25 : ellipse->radius = PETSC_MAX_REAL;
291 25 : ellipse->vscale = 1.0;
292 25 : rg->data = (void*)ellipse;
293 :
294 25 : rg->ops->istrivial = RGIsTrivial_Ellipse;
295 25 : rg->ops->computecontour = RGComputeContour_Ellipse;
296 25 : rg->ops->computebbox = RGComputeBoundingBox_Ellipse;
297 25 : rg->ops->computequadrature = RGComputeQuadrature_Ellipse;
298 25 : rg->ops->checkinside = RGCheckInside_Ellipse;
299 25 : rg->ops->isaxisymmetric = RGIsAxisymmetric_Ellipse;
300 25 : rg->ops->setfromoptions = RGSetFromOptions_Ellipse;
301 25 : rg->ops->view = RGView_Ellipse;
302 25 : rg->ops->destroy = RGDestroy_Ellipse;
303 25 : PetscCall(PetscObjectComposeFunction((PetscObject)rg,"RGEllipseSetParameters_C",RGEllipseSetParameters_Ellipse));
304 25 : PetscCall(PetscObjectComposeFunction((PetscObject)rg,"RGEllipseGetParameters_C",RGEllipseGetParameters_Ellipse));
305 25 : PetscFunctionReturn(PETSC_SUCCESS);
306 : }
|
