Actual source code: rgellipse.c
slepc-main 2024-11-09
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: Region enclosed in an ellipse (aligned with the coordinate axes)
12: */
14: #include <slepc/private/rgimpl.h>
15: #include <petscdraw.h>
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;
23: static PetscErrorCode RGEllipseSetParameters_Ellipse(RG rg,PetscScalar center,PetscReal radius,PetscReal vscale)
24: {
25: RG_ELLIPSE *ctx = (RG_ELLIPSE*)rg->data;
27: PetscFunctionBegin;
28: ctx->center = center;
29: if (radius == (PetscReal)PETSC_DETERMINE) {
30: ctx->radius = 1.0;
31: } else if (radius != (PetscReal)PETSC_CURRENT) {
32: PetscCheck(radius>0.0,PetscObjectComm((PetscObject)rg),PETSC_ERR_ARG_OUTOFRANGE,"The radius argument must be > 0.0");
33: ctx->radius = radius;
34: }
35: if (vscale == (PetscReal)PETSC_DETERMINE) {
36: ctx->vscale = 1.0;
37: } else if (vscale != (PetscReal)PETSC_CURRENT) {
38: PetscCheck(vscale>0.0,PetscObjectComm((PetscObject)rg),PETSC_ERR_ARG_OUTOFRANGE,"The vscale argument must be > 0.0");
39: ctx->vscale = vscale;
40: }
41: PetscFunctionReturn(PETSC_SUCCESS);
42: }
44: /*@
45: RGEllipseSetParameters - Sets the parameters defining the ellipse region.
47: Logically Collective
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
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
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
65: When PETSc is built with real scalars, the center is restricted to a real value.
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.
70: Level: advanced
72: .seealso: RGEllipseGetParameters()
73: @*/
74: PetscErrorCode RGEllipseSetParameters(RG rg,PetscScalar center,PetscReal radius,PetscReal vscale)
75: {
76: PetscFunctionBegin;
81: PetscTryMethod(rg,"RGEllipseSetParameters_C",(RG,PetscScalar,PetscReal,PetscReal),(rg,center,radius,vscale));
82: PetscFunctionReturn(PETSC_SUCCESS);
83: }
85: static PetscErrorCode RGEllipseGetParameters_Ellipse(RG rg,PetscScalar *center,PetscReal *radius,PetscReal *vscale)
86: {
87: RG_ELLIPSE *ctx = (RG_ELLIPSE*)rg->data;
89: PetscFunctionBegin;
90: if (center) *center = ctx->center;
91: if (radius) *radius = ctx->radius;
92: if (vscale) *vscale = ctx->vscale;
93: PetscFunctionReturn(PETSC_SUCCESS);
94: }
96: /*@
97: RGEllipseGetParameters - Gets the parameters that define the ellipse region.
99: Not Collective
101: Input Parameter:
102: . rg - the region context
104: Output Parameters:
105: + center - center of the region
106: . radius - radius of the region
107: - vscale - vertical scale of the region
109: Level: advanced
111: .seealso: RGEllipseSetParameters()
112: @*/
113: PetscErrorCode RGEllipseGetParameters(RG rg,PetscScalar *center,PetscReal *radius,PetscReal *vscale)
114: {
115: PetscFunctionBegin;
117: PetscUseMethod(rg,"RGEllipseGetParameters_C",(RG,PetscScalar*,PetscReal*,PetscReal*),(rg,center,radius,vscale));
118: PetscFunctionReturn(PETSC_SUCCESS);
119: }
121: static PetscErrorCode RGView_Ellipse(RG rg,PetscViewer viewer)
122: {
123: RG_ELLIPSE *ctx = (RG_ELLIPSE*)rg->data;
124: PetscBool isdraw,isascii;
125: int winw,winh;
126: PetscDraw draw;
127: PetscDrawAxis axis;
128: PetscReal cx,cy,r,ab,cd,lx,ly,w,scale=1.2;
129: char str[50];
131: PetscFunctionBegin;
132: PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw));
133: PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
134: if (isascii) {
135: PetscCall(SlepcSNPrintfScalar(str,sizeof(str),ctx->center,PETSC_FALSE));
136: PetscCall(PetscViewerASCIIPrintf(viewer," center: %s, radius: %g, vscale: %g\n",str,RGShowReal(ctx->radius),RGShowReal(ctx->vscale)));
137: } else if (isdraw) {
138: PetscCall(PetscViewerDrawGetDraw(viewer,0,&draw));
139: PetscCall(PetscDrawCheckResizedWindow(draw));
140: PetscCall(PetscDrawGetWindowSize(draw,&winw,&winh));
141: winw = PetscMax(winw,1); winh = PetscMax(winh,1);
142: PetscCall(PetscDrawClear(draw));
143: PetscCall(PetscDrawSetTitle(draw,"Ellipse region"));
144: PetscCall(PetscDrawAxisCreate(draw,&axis));
145: cx = PetscRealPart(ctx->center)*rg->sfactor;
146: cy = PetscImaginaryPart(ctx->center)*rg->sfactor;
147: r = ctx->radius*rg->sfactor;
148: lx = 2*r;
149: ly = 2*r*ctx->vscale;
150: ab = cx;
151: cd = cy;
152: w = scale*PetscMax(lx/winw,ly/winh)/2;
153: PetscCall(PetscDrawAxisSetLimits(axis,ab-w*winw,ab+w*winw,cd-w*winh,cd+w*winh));
154: PetscCall(PetscDrawAxisDraw(axis));
155: PetscCall(PetscDrawAxisDestroy(&axis));
156: PetscCall(PetscDrawEllipse(draw,cx,cy,2*r,2*r*ctx->vscale,PETSC_DRAW_RED));
157: PetscCall(PetscDrawFlush(draw));
158: PetscCall(PetscDrawSave(draw));
159: PetscCall(PetscDrawPause(draw));
160: }
161: PetscFunctionReturn(PETSC_SUCCESS);
162: }
164: static PetscErrorCode RGIsTrivial_Ellipse(RG rg,PetscBool *trivial)
165: {
166: RG_ELLIPSE *ctx = (RG_ELLIPSE*)rg->data;
168: PetscFunctionBegin;
169: if (rg->complement) *trivial = PetscNot(ctx->radius);
170: else *trivial = PetscNot(ctx->radius<PETSC_MAX_REAL);
171: PetscFunctionReturn(PETSC_SUCCESS);
172: }
174: static PetscErrorCode RGComputeContour_Ellipse(RG rg,PetscInt n,PetscScalar *cr,PetscScalar *ci)
175: {
176: RG_ELLIPSE *ctx = (RG_ELLIPSE*)rg->data;
177: PetscReal theta;
178: PetscInt i;
180: PetscFunctionBegin;
181: for (i=0;i<n;i++) {
182: 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: 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: PetscFunctionReturn(PETSC_SUCCESS);
191: }
193: static PetscErrorCode RGComputeBoundingBox_Ellipse(RG rg,PetscReal *a,PetscReal *b,PetscReal *c,PetscReal *d)
194: {
195: RG_ELLIPSE *ctx = (RG_ELLIPSE*)rg->data;
197: PetscFunctionBegin;
198: if (a) *a = PetscRealPart(ctx->center) - ctx->radius;
199: if (b) *b = PetscRealPart(ctx->center) + ctx->radius;
200: if (c) *c = PetscImaginaryPart(ctx->center) - ctx->radius*ctx->vscale;
201: if (d) *d = PetscImaginaryPart(ctx->center) + ctx->radius*ctx->vscale;
202: PetscFunctionReturn(PETSC_SUCCESS);
203: }
205: static PetscErrorCode RGComputeQuadrature_Ellipse(RG rg,RGQuadRule quad,PetscInt n,PetscScalar *z,PetscScalar *zn,PetscScalar *w)
206: {
207: RG_ELLIPSE *ctx = (RG_ELLIPSE*)rg->data;
208: PetscReal theta;
209: PetscInt i;
211: PetscFunctionBegin;
212: 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: 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: PetscFunctionReturn(PETSC_SUCCESS);
225: }
227: static PetscErrorCode RGCheckInside_Ellipse(RG rg,PetscReal px,PetscReal py,PetscInt *inside)
228: {
229: RG_ELLIPSE *ctx = (RG_ELLIPSE*)rg->data;
230: PetscReal dx,dy,r;
232: 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: dx = (px-ctx->center)/ctx->radius;
238: dy = py/ctx->radius;
239: #endif
240: r = 1.0-dx*dx-(dy*dy)/(ctx->vscale*ctx->vscale);
241: *inside = PetscSign(r);
242: PetscFunctionReturn(PETSC_SUCCESS);
243: }
245: static PetscErrorCode RGIsAxisymmetric_Ellipse(RG rg,PetscBool vertical,PetscBool *symm)
246: {
247: RG_ELLIPSE *ctx = (RG_ELLIPSE*)rg->data;
249: PetscFunctionBegin;
250: if (vertical) *symm = (PetscRealPart(ctx->center) == 0.0)? PETSC_TRUE: PETSC_FALSE;
251: else *symm = (PetscImaginaryPart(ctx->center) == 0.0)? PETSC_TRUE: PETSC_FALSE;
252: PetscFunctionReturn(PETSC_SUCCESS);
253: }
255: static PetscErrorCode RGSetFromOptions_Ellipse(RG rg,PetscOptionItems *PetscOptionsObject)
256: {
257: PetscScalar s;
258: PetscReal r1,r2;
259: PetscBool flg1,flg2,flg3;
261: PetscFunctionBegin;
262: PetscOptionsHeadBegin(PetscOptionsObject,"RG Ellipse Options");
264: PetscCall(RGEllipseGetParameters(rg,&s,&r1,&r2));
265: PetscCall(PetscOptionsScalar("-rg_ellipse_center","Center of ellipse","RGEllipseSetParameters",s,&s,&flg1));
266: PetscCall(PetscOptionsReal("-rg_ellipse_radius","Radius of ellipse","RGEllipseSetParameters",r1,&r1,&flg2));
267: PetscCall(PetscOptionsReal("-rg_ellipse_vscale","Vertical scale of ellipse","RGEllipseSetParameters",r2,&r2,&flg3));
268: if (flg1 || flg2 || flg3) PetscCall(RGEllipseSetParameters(rg,s,r1,r2));
270: PetscOptionsHeadEnd();
271: PetscFunctionReturn(PETSC_SUCCESS);
272: }
274: static PetscErrorCode RGDestroy_Ellipse(RG rg)
275: {
276: PetscFunctionBegin;
277: PetscCall(PetscFree(rg->data));
278: PetscCall(PetscObjectComposeFunction((PetscObject)rg,"RGEllipseSetParameters_C",NULL));
279: PetscCall(PetscObjectComposeFunction((PetscObject)rg,"RGEllipseGetParameters_C",NULL));
280: PetscFunctionReturn(PETSC_SUCCESS);
281: }
283: SLEPC_EXTERN PetscErrorCode RGCreate_Ellipse(RG rg)
284: {
285: RG_ELLIPSE *ellipse;
287: PetscFunctionBegin;
288: PetscCall(PetscNew(&ellipse));
289: ellipse->center = 0.0;
290: ellipse->radius = PETSC_MAX_REAL;
291: ellipse->vscale = 1.0;
292: rg->data = (void*)ellipse;
294: rg->ops->istrivial = RGIsTrivial_Ellipse;
295: rg->ops->computecontour = RGComputeContour_Ellipse;
296: rg->ops->computebbox = RGComputeBoundingBox_Ellipse;
297: rg->ops->computequadrature = RGComputeQuadrature_Ellipse;
298: rg->ops->checkinside = RGCheckInside_Ellipse;
299: rg->ops->isaxisymmetric = RGIsAxisymmetric_Ellipse;
300: rg->ops->setfromoptions = RGSetFromOptions_Ellipse;
301: rg->ops->view = RGView_Ellipse;
302: rg->ops->destroy = RGDestroy_Ellipse;
303: PetscCall(PetscObjectComposeFunction((PetscObject)rg,"RGEllipseSetParameters_C",RGEllipseSetParameters_Ellipse));
304: PetscCall(PetscObjectComposeFunction((PetscObject)rg,"RGEllipseGetParameters_C",RGEllipseGetParameters_Ellipse));
305: PetscFunctionReturn(PETSC_SUCCESS);
306: }