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 : Polygonal region defined by a set of vertices
12 : */
13 :
14 : #include <slepc/private/rgimpl.h> /*I "slepcrg.h" I*/
15 : #include <petscdraw.h>
16 :
17 : #define VERTMAX 30
18 :
19 : typedef struct {
20 : PetscInt n; /* number of vertices */
21 : PetscScalar *vr,*vi; /* array of vertices (vi not used in complex scalars) */
22 : } RG_POLYGON;
23 :
24 : static PetscErrorCode RGComputeBoundingBox_Polygon(RG,PetscReal*,PetscReal*,PetscReal*,PetscReal*);
25 :
26 : #if !defined(PETSC_USE_COMPLEX)
27 : static PetscBool CheckSymmetry(PetscInt n,PetscScalar *vr,PetscScalar *vi)
28 : {
29 : PetscInt i,j,k;
30 : /* find change of sign in imaginary part */
31 : j = vi[0]!=0.0? 0: 1;
32 : for (k=j+1;k<n;k++) {
33 : if (vi[k]!=0.0) {
34 : if (vi[k]*vi[j]<0.0) break;
35 : j++;
36 : }
37 : }
38 : if (k==n) return (j==1)? PETSC_TRUE: PETSC_FALSE;
39 : /* check pairing vertices */
40 : for (i=0;i<n/2;i++) {
41 : if (vr[k]!=vr[j] || vi[k]!=-vi[j]) return PETSC_FALSE;
42 : k = (k+1)%n;
43 : j = (j-1+n)%n;
44 : }
45 : return PETSC_TRUE;
46 : }
47 : #endif
48 :
49 3 : static PetscErrorCode RGPolygonSetVertices_Polygon(RG rg,PetscInt n,PetscScalar *vr,PetscScalar *vi)
50 : {
51 3 : PetscInt i;
52 3 : RG_POLYGON *ctx = (RG_POLYGON*)rg->data;
53 :
54 3 : PetscFunctionBegin;
55 3 : PetscCheck(n>2,PetscObjectComm((PetscObject)rg),PETSC_ERR_ARG_OUTOFRANGE,"At least 3 vertices required, you provided %" PetscInt_FMT,n);
56 3 : PetscCheck(n<=VERTMAX,PetscObjectComm((PetscObject)rg),PETSC_ERR_ARG_OUTOFRANGE,"Too many points, maximum allowed is %d",VERTMAX);
57 : #if !defined(PETSC_USE_COMPLEX)
58 : PetscCheck(CheckSymmetry(n,vr,vi),PetscObjectComm((PetscObject)rg),PETSC_ERR_ARG_WRONG,"In real scalars the region must be symmetric wrt real axis");
59 : #endif
60 3 : if (ctx->n) {
61 0 : PetscCall(PetscFree(ctx->vr));
62 : #if !defined(PETSC_USE_COMPLEX)
63 : PetscCall(PetscFree(ctx->vi));
64 : #endif
65 : }
66 3 : ctx->n = n;
67 3 : PetscCall(PetscMalloc1(n,&ctx->vr));
68 : #if !defined(PETSC_USE_COMPLEX)
69 : PetscCall(PetscMalloc1(n,&ctx->vi));
70 : #endif
71 21 : for (i=0;i<n;i++) {
72 18 : ctx->vr[i] = vr[i];
73 : #if !defined(PETSC_USE_COMPLEX)
74 : ctx->vi[i] = vi[i];
75 : #endif
76 : }
77 3 : PetscFunctionReturn(PETSC_SUCCESS);
78 : }
79 :
80 : /*@
81 : RGPolygonSetVertices - Sets the vertices that define the polygon region.
82 :
83 : Logically Collective
84 :
85 : Input Parameters:
86 : + rg - the region context
87 : . n - number of vertices
88 : . vr - array of vertices
89 : - vi - array of vertices (imaginary part)
90 :
91 : Options Database Keys:
92 : + -rg_polygon_vertices - Sets the vertices
93 : - -rg_polygon_verticesi - Sets the vertices (imaginary part)
94 :
95 : Notes:
96 : In the case of complex scalars, only argument vr is used, containing
97 : the complex vertices; the list of vertices can be provided in the
98 : command line with a comma-separated list of complex values
99 : [+/-][realnumber][+/-]realnumberi with no spaces.
100 :
101 : When PETSc is built with real scalars, the real and imaginary parts of
102 : the vertices must be provided in two separate arrays (or two lists in
103 : the command line). In this case, the region must be symmetric with
104 : respect to the real axis.
105 :
106 : Level: advanced
107 :
108 : .seealso: RGPolygonGetVertices()
109 : @*/
110 3 : PetscErrorCode RGPolygonSetVertices(RG rg,PetscInt n,PetscScalar vr[],PetscScalar vi[])
111 : {
112 3 : PetscFunctionBegin;
113 3 : PetscValidHeaderSpecific(rg,RG_CLASSID,1);
114 12 : PetscValidLogicalCollectiveInt(rg,n,2);
115 3 : PetscAssertPointer(vr,3);
116 : #if !defined(PETSC_USE_COMPLEX)
117 : PetscAssertPointer(vi,4);
118 : #endif
119 3 : PetscTryMethod(rg,"RGPolygonSetVertices_C",(RG,PetscInt,PetscScalar*,PetscScalar*),(rg,n,vr,vi));
120 3 : PetscFunctionReturn(PETSC_SUCCESS);
121 : }
122 :
123 2 : static PetscErrorCode RGPolygonGetVertices_Polygon(RG rg,PetscInt *n,PetscScalar **vr,PetscScalar **vi)
124 : {
125 2 : RG_POLYGON *ctx = (RG_POLYGON*)rg->data;
126 2 : PetscInt i;
127 :
128 2 : PetscFunctionBegin;
129 2 : if (n) *n = ctx->n;
130 2 : if (vr) {
131 2 : if (!ctx->n) *vr = NULL;
132 : else {
133 2 : PetscCall(PetscMalloc1(ctx->n,vr));
134 16 : for (i=0;i<ctx->n;i++) (*vr)[i] = ctx->vr[i];
135 : }
136 : }
137 : #if !defined(PETSC_USE_COMPLEX)
138 : if (vi) {
139 : if (!ctx->n) *vi = NULL;
140 : else {
141 : PetscCall(PetscMalloc1(ctx->n,vi));
142 : for (i=0;i<ctx->n;i++) (*vi)[i] = ctx->vi[i];
143 : }
144 : }
145 : #endif
146 2 : PetscFunctionReturn(PETSC_SUCCESS);
147 : }
148 :
149 : /*@C
150 : RGPolygonGetVertices - Gets the vertices that define the polygon region.
151 :
152 : Not Collective
153 :
154 : Input Parameter:
155 : . rg - the region context
156 :
157 : Output Parameters:
158 : + n - number of vertices
159 : . vr - array of vertices
160 : - vi - array of vertices (imaginary part)
161 :
162 : Notes:
163 : The values passed by user with RGPolygonSetVertices() are returned (or null
164 : pointers otherwise).
165 : The returned arrays should be freed by the user when no longer needed.
166 :
167 : Level: advanced
168 :
169 : .seealso: RGPolygonSetVertices()
170 : @*/
171 2 : PetscErrorCode RGPolygonGetVertices(RG rg,PetscInt *n,PetscScalar **vr,PetscScalar **vi)
172 : {
173 2 : PetscFunctionBegin;
174 2 : PetscValidHeaderSpecific(rg,RG_CLASSID,1);
175 2 : PetscUseMethod(rg,"RGPolygonGetVertices_C",(RG,PetscInt*,PetscScalar**,PetscScalar**),(rg,n,vr,vi));
176 2 : PetscFunctionReturn(PETSC_SUCCESS);
177 : }
178 :
179 3 : static PetscErrorCode RGView_Polygon(RG rg,PetscViewer viewer)
180 : {
181 3 : RG_POLYGON *ctx = (RG_POLYGON*)rg->data;
182 3 : PetscBool isdraw,isascii;
183 3 : int winw,winh;
184 3 : PetscDraw draw;
185 3 : PetscDrawAxis axis;
186 3 : PetscReal a,b,c,d,ab,cd,lx,ly,w,x0,y0,x1,y1,scale=1.2;
187 3 : PetscInt i;
188 3 : char str[50];
189 :
190 3 : PetscFunctionBegin;
191 3 : PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw));
192 3 : PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
193 3 : if (isascii) {
194 2 : PetscCall(PetscViewerASCIIPrintf(viewer," vertices: "));
195 2 : PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_FALSE));
196 16 : for (i=0;i<ctx->n;i++) {
197 : #if defined(PETSC_USE_COMPLEX)
198 14 : PetscCall(SlepcSNPrintfScalar(str,sizeof(str),ctx->vr[i],PETSC_FALSE));
199 : #else
200 : if (ctx->vi[i]!=0.0) PetscCall(PetscSNPrintf(str,sizeof(str),"%g%+gi",(double)ctx->vr[i],(double)ctx->vi[i]));
201 : else PetscCall(PetscSNPrintf(str,sizeof(str),"%g",(double)ctx->vr[i]));
202 : #endif
203 16 : PetscCall(PetscViewerASCIIPrintf(viewer,"%s%s",str,(i<ctx->n-1)?", ":""));
204 : }
205 2 : PetscCall(PetscViewerASCIIPrintf(viewer,"\n"));
206 2 : PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_TRUE));
207 1 : } else if (isdraw) {
208 1 : PetscCall(PetscViewerDrawGetDraw(viewer,0,&draw));
209 1 : PetscCall(PetscDrawCheckResizedWindow(draw));
210 1 : PetscCall(PetscDrawGetWindowSize(draw,&winw,&winh));
211 1 : winw = PetscMax(winw,1); winh = PetscMax(winh,1);
212 1 : PetscCall(PetscDrawClear(draw));
213 1 : PetscCall(PetscDrawSetTitle(draw,"Polygonal region"));
214 1 : PetscCall(PetscDrawAxisCreate(draw,&axis));
215 1 : PetscCall(RGComputeBoundingBox_Polygon(rg,&a,&b,&c,&d));
216 1 : a *= rg->sfactor;
217 1 : b *= rg->sfactor;
218 1 : c *= rg->sfactor;
219 1 : d *= rg->sfactor;
220 1 : lx = b-a;
221 1 : ly = d-c;
222 1 : ab = (a+b)/2;
223 1 : cd = (c+d)/2;
224 1 : w = scale*PetscMax(lx/winw,ly/winh)/2;
225 1 : PetscCall(PetscDrawAxisSetLimits(axis,ab-w*winw,ab+w*winw,cd-w*winh,cd+w*winh));
226 1 : PetscCall(PetscDrawAxisDraw(axis));
227 1 : PetscCall(PetscDrawAxisDestroy(&axis));
228 8 : for (i=0;i<ctx->n;i++) {
229 : #if defined(PETSC_USE_COMPLEX)
230 7 : x0 = PetscRealPart(ctx->vr[i]); y0 = PetscImaginaryPart(ctx->vr[i]);
231 7 : if (i<ctx->n-1) {
232 6 : x1 = PetscRealPart(ctx->vr[i+1]); y1 = PetscImaginaryPart(ctx->vr[i+1]);
233 : } else {
234 1 : x1 = PetscRealPart(ctx->vr[0]); y1 = PetscImaginaryPart(ctx->vr[0]);
235 : }
236 : #else
237 : x0 = ctx->vr[i]; y0 = ctx->vi[i];
238 : if (i<ctx->n-1) {
239 : x1 = ctx->vr[i+1]; y1 = ctx->vi[i+1];
240 : } else {
241 : x1 = ctx->vr[0]; y1 = ctx->vi[0];
242 : }
243 : #endif
244 7 : PetscCall(PetscDrawLine(draw,x0*rg->sfactor,y0*rg->sfactor,x1*rg->sfactor,y1*rg->sfactor,PETSC_DRAW_MAGENTA));
245 : }
246 1 : PetscCall(PetscDrawFlush(draw));
247 1 : PetscCall(PetscDrawSave(draw));
248 1 : PetscCall(PetscDrawPause(draw));
249 : }
250 3 : PetscFunctionReturn(PETSC_SUCCESS);
251 : }
252 :
253 10 : static PetscErrorCode RGIsTrivial_Polygon(RG rg,PetscBool *trivial)
254 : {
255 10 : RG_POLYGON *ctx = (RG_POLYGON*)rg->data;
256 :
257 10 : PetscFunctionBegin;
258 10 : *trivial = PetscNot(ctx->n);
259 10 : PetscFunctionReturn(PETSC_SUCCESS);
260 : }
261 :
262 3 : static PetscErrorCode RGComputeContour_Polygon(RG rg,PetscInt n,PetscScalar *ucr,PetscScalar *uci)
263 : {
264 3 : RG_POLYGON *ctx = (RG_POLYGON*)rg->data;
265 3 : PetscReal length,h,d,rem=0.0;
266 3 : PetscInt k=1,idx=ctx->n-1,i;
267 3 : PetscBool ini=PETSC_FALSE;
268 3 : PetscScalar incr,*cr=ucr,*ci=uci;
269 : #if !defined(PETSC_USE_COMPLEX)
270 : PetscScalar inci;
271 : #endif
272 :
273 3 : PetscFunctionBegin;
274 3 : PetscCheck(ctx->n,PetscObjectComm((PetscObject)rg),PETSC_ERR_ARG_WRONGSTATE,"No vertices have been set yet");
275 3 : length = SlepcAbsEigenvalue(ctx->vr[0]-ctx->vr[ctx->n-1],ctx->vi[0]-ctx->vi[ctx->n-1]);
276 18 : for (i=0;i<ctx->n-1;i++) length += SlepcAbsEigenvalue(ctx->vr[i]-ctx->vr[i+1],ctx->vi[i]-ctx->vi[i+1]);
277 3 : h = length/n;
278 3 : if (!ucr) PetscCall(PetscMalloc1(n,&cr));
279 3 : if (!uci) PetscCall(PetscMalloc1(n,&ci));
280 3 : cr[0] = ctx->vr[0];
281 : #if !defined(PETSC_USE_COMPLEX)
282 : ci[0] = ctx->vi[0];
283 : #endif
284 3 : incr = ctx->vr[ctx->n-1]-ctx->vr[0];
285 : #if !defined(PETSC_USE_COMPLEX)
286 : inci = ctx->vi[ctx->n-1]-ctx->vi[0];
287 : #endif
288 3 : d = SlepcAbsEigenvalue(incr,inci);
289 3 : incr /= d;
290 : #if !defined(PETSC_USE_COMPLEX)
291 : inci /= d;
292 : #endif
293 135 : while (k<n) {
294 132 : if (ini) {
295 15 : incr = ctx->vr[idx]-ctx->vr[idx+1];
296 : #if !defined(PETSC_USE_COMPLEX)
297 : inci = ctx->vi[idx]-ctx->vi[idx+1];
298 : #endif
299 15 : d = SlepcAbsEigenvalue(incr,inci);
300 15 : incr /= d;
301 : #if !defined(PETSC_USE_COMPLEX)
302 : inci /= d;
303 : #endif
304 15 : if (rem+d>h) {
305 15 : cr[k] = ctx->vr[idx+1]+incr*(h-rem);
306 : #if !defined(PETSC_USE_COMPLEX)
307 : ci[k] = ctx->vi[idx+1]+inci*(h-rem);
308 : #endif
309 15 : k++;
310 15 : ini = PETSC_FALSE;
311 0 : } else {rem += d; idx--;}
312 : } else {
313 : #if !defined(PETSC_USE_COMPLEX)
314 : rem = SlepcAbsEigenvalue(ctx->vr[idx]-cr[k-1],ctx->vi[idx]-ci[k-1]);
315 : #else
316 117 : rem = PetscAbsScalar(ctx->vr[idx]-cr[k-1]);
317 : #endif
318 117 : if (rem>h) {
319 102 : cr[k] = cr[k-1]+incr*h;
320 : #if !defined(PETSC_USE_COMPLEX)
321 : ci[k] = ci[k-1]+inci*h;
322 : #endif
323 102 : k++;
324 15 : } else {ini = PETSC_TRUE; idx--;}
325 : }
326 : }
327 3 : if (!ucr) PetscCall(PetscFree(cr));
328 3 : if (!uci) PetscCall(PetscFree(ci));
329 3 : PetscFunctionReturn(PETSC_SUCCESS);
330 : }
331 :
332 3 : static PetscErrorCode RGComputeBoundingBox_Polygon(RG rg,PetscReal *a,PetscReal *b,PetscReal *c,PetscReal *d)
333 : {
334 3 : RG_POLYGON *ctx = (RG_POLYGON*)rg->data;
335 3 : PetscInt i;
336 :
337 3 : PetscFunctionBegin;
338 3 : if (a) *a = PETSC_MAX_REAL;
339 3 : if (b) *b = -PETSC_MAX_REAL;
340 3 : if (c) *c = PETSC_MAX_REAL;
341 3 : if (d) *d = -PETSC_MAX_REAL;
342 24 : for (i=0;i<ctx->n;i++) {
343 : #if defined(PETSC_USE_COMPLEX)
344 24 : if (a) *a = PetscMin(*a,PetscRealPart(ctx->vr[i]));
345 24 : if (b) *b = PetscMax(*b,PetscRealPart(ctx->vr[i]));
346 30 : if (c) *c = PetscMin(*c,PetscImaginaryPart(ctx->vr[i]));
347 33 : if (d) *d = PetscMax(*d,PetscImaginaryPart(ctx->vr[i]));
348 : #else
349 : if (a) *a = PetscMin(*a,ctx->vr[i]);
350 : if (b) *b = PetscMax(*b,ctx->vr[i]);
351 : if (c) *c = PetscMin(*c,ctx->vi[i]);
352 : if (d) *d = PetscMax(*d,ctx->vi[i]);
353 : #endif
354 : }
355 3 : PetscFunctionReturn(PETSC_SUCCESS);
356 : }
357 :
358 5946 : static PetscErrorCode RGCheckInside_Polygon(RG rg,PetscReal px,PetscReal py,PetscInt *inout)
359 : {
360 5946 : RG_POLYGON *ctx = (RG_POLYGON*)rg->data;
361 5946 : PetscReal val,x[VERTMAX],y[VERTMAX];
362 5946 : PetscBool mx,my,nx,ny;
363 5946 : PetscInt i,j;
364 :
365 5946 : PetscFunctionBegin;
366 29736 : for (i=0;i<ctx->n;i++) {
367 : #if defined(PETSC_USE_COMPLEX)
368 23790 : x[i] = PetscRealPart(ctx->vr[i])-px;
369 23790 : y[i] = PetscImaginaryPart(ctx->vr[i])-py;
370 : #else
371 : x[i] = ctx->vr[i]-px;
372 : y[i] = ctx->vi[i]-py;
373 : #endif
374 : }
375 5946 : *inout = -1;
376 29736 : for (i=0;i<ctx->n;i++) {
377 23790 : j = (i+1)%ctx->n;
378 23790 : mx = PetscNot(x[i]<0.0);
379 23790 : nx = PetscNot(x[j]<0.0);
380 23790 : my = PetscNot(y[i]<0.0);
381 23790 : ny = PetscNot(y[j]<0.0);
382 23790 : if (!((my||ny) && (mx||nx)) || (mx&&nx)) continue;
383 2267 : if (((my && ny && (mx||nx)) && (!(mx&&nx)))) {
384 1957 : *inout = -*inout;
385 1957 : continue;
386 : }
387 310 : val = (y[i]*x[j]-x[i]*y[j])/(x[j]-x[i]);
388 310 : if (PetscAbs(val)<10*PETSC_MACHINE_EPSILON) {
389 0 : *inout = 0;
390 0 : PetscFunctionReturn(PETSC_SUCCESS);
391 310 : } else if (val>0.0) *inout = -*inout;
392 : }
393 5946 : PetscFunctionReturn(PETSC_SUCCESS);
394 : }
395 :
396 3 : static PetscErrorCode RGSetFromOptions_Polygon(RG rg,PetscOptionItems *PetscOptionsObject)
397 : {
398 3 : PetscScalar array[VERTMAX];
399 3 : PetscInt i,k;
400 3 : PetscBool flg,flgi=PETSC_FALSE;
401 : #if !defined(PETSC_USE_COMPLEX)
402 : PetscScalar arrayi[VERTMAX];
403 : PetscInt ki;
404 : #else
405 3 : PetscScalar *arrayi=NULL;
406 : #endif
407 :
408 3 : PetscFunctionBegin;
409 3 : PetscOptionsHeadBegin(PetscOptionsObject,"RG Polygon Options");
410 :
411 3 : k = VERTMAX;
412 93 : for (i=0;i<k;i++) array[i] = 0;
413 3 : PetscCall(PetscOptionsScalarArray("-rg_polygon_vertices","Vertices of polygon","RGPolygonSetVertices",array,&k,&flg));
414 : #if !defined(PETSC_USE_COMPLEX)
415 : ki = VERTMAX;
416 : for (i=0;i<ki;i++) arrayi[i] = 0;
417 : PetscCall(PetscOptionsScalarArray("-rg_polygon_verticesi","Vertices of polygon (imaginary part)","RGPolygonSetVertices",arrayi,&ki,&flgi));
418 : PetscCheck(ki==k,PetscObjectComm((PetscObject)rg),PETSC_ERR_ARG_SIZ,"The number of real %" PetscInt_FMT " and imaginary %" PetscInt_FMT " parts do not match",k,ki);
419 : #endif
420 3 : if (flg || flgi) PetscCall(RGPolygonSetVertices(rg,k,array,arrayi));
421 :
422 3 : PetscOptionsHeadEnd();
423 3 : PetscFunctionReturn(PETSC_SUCCESS);
424 : }
425 :
426 3 : static PetscErrorCode RGDestroy_Polygon(RG rg)
427 : {
428 3 : RG_POLYGON *ctx = (RG_POLYGON*)rg->data;
429 :
430 3 : PetscFunctionBegin;
431 3 : if (ctx->n) {
432 3 : PetscCall(PetscFree(ctx->vr));
433 : #if !defined(PETSC_USE_COMPLEX)
434 : PetscCall(PetscFree(ctx->vi));
435 : #endif
436 : }
437 3 : PetscCall(PetscFree(rg->data));
438 3 : PetscCall(PetscObjectComposeFunction((PetscObject)rg,"RGPolygonSetVertices_C",NULL));
439 3 : PetscCall(PetscObjectComposeFunction((PetscObject)rg,"RGPolygonGetVertices_C",NULL));
440 3 : PetscFunctionReturn(PETSC_SUCCESS);
441 : }
442 :
443 3 : SLEPC_EXTERN PetscErrorCode RGCreate_Polygon(RG rg)
444 : {
445 3 : RG_POLYGON *polygon;
446 :
447 3 : PetscFunctionBegin;
448 3 : PetscCall(PetscNew(&polygon));
449 3 : rg->data = (void*)polygon;
450 :
451 3 : rg->ops->istrivial = RGIsTrivial_Polygon;
452 3 : rg->ops->computecontour = RGComputeContour_Polygon;
453 3 : rg->ops->computebbox = RGComputeBoundingBox_Polygon;
454 3 : rg->ops->checkinside = RGCheckInside_Polygon;
455 3 : rg->ops->setfromoptions = RGSetFromOptions_Polygon;
456 3 : rg->ops->view = RGView_Polygon;
457 3 : rg->ops->destroy = RGDestroy_Polygon;
458 3 : PetscCall(PetscObjectComposeFunction((PetscObject)rg,"RGPolygonSetVertices_C",RGPolygonSetVertices_Polygon));
459 3 : PetscCall(PetscObjectComposeFunction((PetscObject)rg,"RGPolygonGetVertices_C",RGPolygonGetVertices_Polygon));
460 3 : PetscFunctionReturn(PETSC_SUCCESS);
461 : }
|
