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 : Basic RG routines
12 : */
13 :
14 : #include <slepc/private/rgimpl.h> /*I "slepcrg.h" I*/
15 :
16 : PetscFunctionList RGList = NULL;
17 : PetscBool RGRegisterAllCalled = PETSC_FALSE;
18 : PetscClassId RG_CLASSID = 0;
19 : static PetscBool RGPackageInitialized = PETSC_FALSE;
20 :
21 : /*@C
22 : RGFinalizePackage - This function destroys everything in the Slepc interface
23 : to the RG package. It is called from SlepcFinalize().
24 :
25 : Level: developer
26 :
27 : .seealso: SlepcFinalize()
28 : @*/
29 786 : PetscErrorCode RGFinalizePackage(void)
30 : {
31 786 : PetscFunctionBegin;
32 786 : PetscCall(PetscFunctionListDestroy(&RGList));
33 786 : RGPackageInitialized = PETSC_FALSE;
34 786 : RGRegisterAllCalled = PETSC_FALSE;
35 786 : PetscFunctionReturn(PETSC_SUCCESS);
36 : }
37 :
38 : /*@C
39 : RGInitializePackage - This function initializes everything in the RG package.
40 : It is called from PetscDLLibraryRegister() when using dynamic libraries, and
41 : on the first call to RGCreate() when using static libraries.
42 :
43 : Level: developer
44 :
45 : .seealso: SlepcInitialize()
46 : @*/
47 4065 : PetscErrorCode RGInitializePackage(void)
48 : {
49 4065 : char logList[256];
50 4065 : PetscBool opt,pkg;
51 4065 : PetscClassId classids[1];
52 :
53 4065 : PetscFunctionBegin;
54 4065 : if (RGPackageInitialized) PetscFunctionReturn(PETSC_SUCCESS);
55 786 : RGPackageInitialized = PETSC_TRUE;
56 : /* Register Classes */
57 786 : PetscCall(PetscClassIdRegister("Region",&RG_CLASSID));
58 : /* Register Constructors */
59 786 : PetscCall(RGRegisterAll());
60 : /* Process Info */
61 786 : classids[0] = RG_CLASSID;
62 786 : PetscCall(PetscInfoProcessClass("rg",1,&classids[0]));
63 : /* Process summary exclusions */
64 786 : PetscCall(PetscOptionsGetString(NULL,NULL,"-log_exclude",logList,sizeof(logList),&opt));
65 786 : if (opt) {
66 8 : PetscCall(PetscStrInList("rg",logList,',',&pkg));
67 8 : if (pkg) PetscCall(PetscLogEventDeactivateClass(RG_CLASSID));
68 : }
69 : /* Register package finalizer */
70 786 : PetscCall(PetscRegisterFinalize(RGFinalizePackage));
71 786 : PetscFunctionReturn(PETSC_SUCCESS);
72 : }
73 :
74 : /*@
75 : RGCreate - Creates an RG context.
76 :
77 : Collective
78 :
79 : Input Parameter:
80 : . comm - MPI communicator
81 :
82 : Output Parameter:
83 : . newrg - location to put the RG context
84 :
85 : Level: beginner
86 :
87 : .seealso: RGDestroy(), RG
88 : @*/
89 920 : PetscErrorCode RGCreate(MPI_Comm comm,RG *newrg)
90 : {
91 920 : RG rg;
92 :
93 920 : PetscFunctionBegin;
94 920 : PetscAssertPointer(newrg,2);
95 920 : PetscCall(RGInitializePackage());
96 920 : PetscCall(SlepcHeaderCreate(rg,RG_CLASSID,"RG","Region","RG",comm,RGDestroy,RGView));
97 920 : rg->complement = PETSC_FALSE;
98 920 : rg->sfactor = 1.0;
99 920 : rg->osfactor = 0.0;
100 920 : rg->data = NULL;
101 :
102 920 : *newrg = rg;
103 920 : PetscFunctionReturn(PETSC_SUCCESS);
104 : }
105 :
106 : /*@
107 : RGSetOptionsPrefix - Sets the prefix used for searching for all
108 : RG options in the database.
109 :
110 : Logically Collective
111 :
112 : Input Parameters:
113 : + rg - the region context
114 : - prefix - the prefix string to prepend to all RG option requests
115 :
116 : Notes:
117 : A hyphen (-) must NOT be given at the beginning of the prefix name.
118 : The first character of all runtime options is AUTOMATICALLY the
119 : hyphen.
120 :
121 : Level: advanced
122 :
123 : .seealso: RGAppendOptionsPrefix()
124 : @*/
125 227 : PetscErrorCode RGSetOptionsPrefix(RG rg,const char *prefix)
126 : {
127 227 : PetscFunctionBegin;
128 227 : PetscValidHeaderSpecific(rg,RG_CLASSID,1);
129 227 : PetscCall(PetscObjectSetOptionsPrefix((PetscObject)rg,prefix));
130 227 : PetscFunctionReturn(PETSC_SUCCESS);
131 : }
132 :
133 : /*@
134 : RGAppendOptionsPrefix - Appends to the prefix used for searching for all
135 : RG options in the database.
136 :
137 : Logically Collective
138 :
139 : Input Parameters:
140 : + rg - the region context
141 : - prefix - the prefix string to prepend to all RG option requests
142 :
143 : Notes:
144 : A hyphen (-) must NOT be given at the beginning of the prefix name.
145 : The first character of all runtime options is AUTOMATICALLY the hyphen.
146 :
147 : Level: advanced
148 :
149 : .seealso: RGSetOptionsPrefix()
150 : @*/
151 187 : PetscErrorCode RGAppendOptionsPrefix(RG rg,const char *prefix)
152 : {
153 187 : PetscFunctionBegin;
154 187 : PetscValidHeaderSpecific(rg,RG_CLASSID,1);
155 187 : PetscCall(PetscObjectAppendOptionsPrefix((PetscObject)rg,prefix));
156 187 : PetscFunctionReturn(PETSC_SUCCESS);
157 : }
158 :
159 : /*@
160 : RGGetOptionsPrefix - Gets the prefix used for searching for all
161 : RG options in the database.
162 :
163 : Not Collective
164 :
165 : Input Parameters:
166 : . rg - the region context
167 :
168 : Output Parameters:
169 : . prefix - pointer to the prefix string used is returned
170 :
171 : Note:
172 : On the Fortran side, the user should pass in a string 'prefix' of
173 : sufficient length to hold the prefix.
174 :
175 : Level: advanced
176 :
177 : .seealso: RGSetOptionsPrefix(), RGAppendOptionsPrefix()
178 : @*/
179 0 : PetscErrorCode RGGetOptionsPrefix(RG rg,const char *prefix[])
180 : {
181 0 : PetscFunctionBegin;
182 0 : PetscValidHeaderSpecific(rg,RG_CLASSID,1);
183 0 : PetscAssertPointer(prefix,2);
184 0 : PetscCall(PetscObjectGetOptionsPrefix((PetscObject)rg,prefix));
185 0 : PetscFunctionReturn(PETSC_SUCCESS);
186 : }
187 :
188 : /*@
189 : RGSetType - Selects the type for the RG object.
190 :
191 : Logically Collective
192 :
193 : Input Parameters:
194 : + rg - the region context
195 : - type - a known type
196 :
197 : Level: intermediate
198 :
199 : .seealso: RGGetType()
200 : @*/
201 966 : PetscErrorCode RGSetType(RG rg,RGType type)
202 : {
203 966 : PetscErrorCode (*r)(RG);
204 966 : PetscBool match;
205 :
206 966 : PetscFunctionBegin;
207 966 : PetscValidHeaderSpecific(rg,RG_CLASSID,1);
208 966 : PetscAssertPointer(type,2);
209 :
210 966 : PetscCall(PetscObjectTypeCompare((PetscObject)rg,type,&match));
211 966 : if (match) PetscFunctionReturn(PETSC_SUCCESS);
212 :
213 924 : PetscCall(PetscFunctionListFind(RGList,type,&r));
214 924 : PetscCheck(r,PetscObjectComm((PetscObject)rg),PETSC_ERR_ARG_UNKNOWN_TYPE,"Unable to find requested RG type %s",type);
215 :
216 924 : PetscTryTypeMethod(rg,destroy);
217 924 : PetscCall(PetscMemzero(rg->ops,sizeof(struct _RGOps)));
218 :
219 924 : PetscCall(PetscObjectChangeTypeName((PetscObject)rg,type));
220 924 : PetscCall((*r)(rg));
221 924 : PetscFunctionReturn(PETSC_SUCCESS);
222 : }
223 :
224 : /*@
225 : RGGetType - Gets the RG type name (as a string) from the RG context.
226 :
227 : Not Collective
228 :
229 : Input Parameter:
230 : . rg - the region context
231 :
232 : Output Parameter:
233 : . type - name of the region
234 :
235 : Level: intermediate
236 :
237 : .seealso: RGSetType()
238 : @*/
239 5 : PetscErrorCode RGGetType(RG rg,RGType *type)
240 : {
241 5 : PetscFunctionBegin;
242 5 : PetscValidHeaderSpecific(rg,RG_CLASSID,1);
243 5 : PetscAssertPointer(type,2);
244 5 : *type = ((PetscObject)rg)->type_name;
245 5 : PetscFunctionReturn(PETSC_SUCCESS);
246 : }
247 :
248 : /*@
249 : RGSetFromOptions - Sets RG options from the options database.
250 :
251 : Collective
252 :
253 : Input Parameters:
254 : . rg - the region context
255 :
256 : Notes:
257 : To see all options, run your program with the -help option.
258 :
259 : Level: beginner
260 :
261 : .seealso: RGSetOptionsPrefix()
262 : @*/
263 898 : PetscErrorCode RGSetFromOptions(RG rg)
264 : {
265 898 : char type[256];
266 898 : PetscBool flg;
267 898 : PetscReal sfactor;
268 :
269 898 : PetscFunctionBegin;
270 898 : PetscValidHeaderSpecific(rg,RG_CLASSID,1);
271 898 : PetscCall(RGRegisterAll());
272 2694 : PetscObjectOptionsBegin((PetscObject)rg);
273 967 : PetscCall(PetscOptionsFList("-rg_type","Region type","RGSetType",RGList,(char*)(((PetscObject)rg)->type_name?((PetscObject)rg)->type_name:RGINTERVAL),type,sizeof(type),&flg));
274 898 : if (flg) PetscCall(RGSetType(rg,type));
275 826 : else if (!((PetscObject)rg)->type_name) PetscCall(RGSetType(rg,RGINTERVAL));
276 :
277 898 : PetscCall(PetscOptionsBool("-rg_complement","Whether region is complemented or not","RGSetComplement",rg->complement,&rg->complement,NULL));
278 :
279 898 : PetscCall(PetscOptionsReal("-rg_scale","Scaling factor","RGSetScale",1.0,&sfactor,&flg));
280 898 : if (flg) PetscCall(RGSetScale(rg,sfactor));
281 :
282 898 : PetscTryTypeMethod(rg,setfromoptions,PetscOptionsObject);
283 898 : PetscCall(PetscObjectProcessOptionsHandlers((PetscObject)rg,PetscOptionsObject));
284 898 : PetscOptionsEnd();
285 898 : PetscFunctionReturn(PETSC_SUCCESS);
286 : }
287 :
288 : /*@
289 : RGView - Prints the RG data structure.
290 :
291 : Collective
292 :
293 : Input Parameters:
294 : + rg - the region context
295 : - viewer - optional visualization context
296 :
297 : Note:
298 : The available visualization contexts include
299 : + PETSC_VIEWER_STDOUT_SELF - standard output (default)
300 : - PETSC_VIEWER_STDOUT_WORLD - synchronized standard
301 : output where only the first processor opens
302 : the file. All other processors send their
303 : data to the first processor to print.
304 :
305 : The user can open an alternative visualization context with
306 : PetscViewerASCIIOpen() - output to a specified file.
307 :
308 : Level: beginner
309 :
310 : .seealso: RGCreate()
311 : @*/
312 17 : PetscErrorCode RGView(RG rg,PetscViewer viewer)
313 : {
314 17 : PetscBool isdraw,isascii;
315 :
316 17 : PetscFunctionBegin;
317 17 : PetscValidHeaderSpecific(rg,RG_CLASSID,1);
318 17 : if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)rg),&viewer));
319 17 : PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,2);
320 17 : PetscCheckSameComm(rg,1,viewer,2);
321 17 : PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw));
322 17 : PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
323 17 : if (isascii) {
324 13 : PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)rg,viewer));
325 13 : PetscCall(PetscViewerASCIIPushTab(viewer));
326 13 : PetscTryTypeMethod(rg,view,viewer);
327 13 : PetscCall(PetscViewerASCIIPopTab(viewer));
328 13 : if (rg->complement) PetscCall(PetscViewerASCIIPrintf(viewer," selected region is the complement of the specified one\n"));
329 13 : if (rg->sfactor!=1.0) PetscCall(PetscViewerASCIIPrintf(viewer," scaling factor = %g\n",(double)rg->sfactor));
330 4 : } else if (isdraw) PetscTryTypeMethod(rg,view,viewer);
331 17 : PetscFunctionReturn(PETSC_SUCCESS);
332 : }
333 :
334 : /*@
335 : RGViewFromOptions - View from options
336 :
337 : Collective
338 :
339 : Input Parameters:
340 : + rg - the region context
341 : . obj - optional object
342 : - name - command line option
343 :
344 : Level: intermediate
345 :
346 : .seealso: RGView(), RGCreate()
347 : @*/
348 1182 : PetscErrorCode RGViewFromOptions(RG rg,PetscObject obj,const char name[])
349 : {
350 1182 : PetscFunctionBegin;
351 1182 : PetscValidHeaderSpecific(rg,RG_CLASSID,1);
352 1182 : PetscCall(PetscObjectViewFromOptions((PetscObject)rg,obj,name));
353 1182 : PetscFunctionReturn(PETSC_SUCCESS);
354 : }
355 :
356 : /*@
357 : RGIsTrivial - Whether it is the trivial region (whole complex plane).
358 :
359 : Not Collective
360 :
361 : Input Parameter:
362 : . rg - the region context
363 :
364 : Output Parameter:
365 : . trivial - true if the region is equal to the whole complex plane, e.g.,
366 : an interval region with all four endpoints unbounded or an
367 : ellipse with infinite radius.
368 :
369 : Level: beginner
370 :
371 : .seealso: RGCheckInside()
372 : @*/
373 9243 : PetscErrorCode RGIsTrivial(RG rg,PetscBool *trivial)
374 : {
375 9243 : PetscFunctionBegin;
376 9243 : PetscValidHeaderSpecific(rg,RG_CLASSID,1);
377 9243 : PetscValidType(rg,1);
378 9243 : PetscAssertPointer(trivial,2);
379 9243 : *trivial = PETSC_FALSE;
380 9243 : PetscTryTypeMethod(rg,istrivial,trivial);
381 9243 : PetscFunctionReturn(PETSC_SUCCESS);
382 : }
383 :
384 : /*@
385 : RGCheckInside - Determines if a set of given points are inside the region or not.
386 :
387 : Not Collective
388 :
389 : Input Parameters:
390 : + rg - the region context
391 : . n - number of points to check
392 : . ar - array of real parts
393 : - ai - array of imaginary parts
394 :
395 : Output Parameter:
396 : . inside - array of results (1=inside, 0=on the contour, -1=outside)
397 :
398 : Note:
399 : The point a is expressed as a couple of PetscScalar variables ar,ai.
400 : If built with complex scalars, the point is supposed to be stored in ar,
401 : otherwise ar,ai contain the real and imaginary parts, respectively.
402 :
403 : If a scaling factor was set, the points are scaled before checking.
404 :
405 : Level: intermediate
406 :
407 : .seealso: RGSetScale(), RGSetComplement()
408 : @*/
409 66178 : PetscErrorCode RGCheckInside(RG rg,PetscInt n,PetscScalar *ar,PetscScalar *ai,PetscInt *inside)
410 : {
411 66178 : PetscReal px,py;
412 66178 : PetscInt i;
413 :
414 66178 : PetscFunctionBegin;
415 66178 : PetscValidHeaderSpecific(rg,RG_CLASSID,1);
416 66178 : PetscValidType(rg,1);
417 66178 : PetscAssertPointer(ar,3);
418 : #if !defined(PETSC_USE_COMPLEX)
419 66178 : PetscAssertPointer(ai,4);
420 : #endif
421 66178 : PetscAssertPointer(inside,5);
422 :
423 198295 : for (i=0;i<n;i++) {
424 : #if defined(PETSC_USE_COMPLEX)
425 : px = PetscRealPart(ar[i]);
426 : py = PetscImaginaryPart(ar[i]);
427 : #else
428 132117 : px = ar[i];
429 132117 : py = ai[i];
430 : #endif
431 132117 : if (PetscUnlikely(rg->sfactor != 1.0)) {
432 2 : px /= rg->sfactor;
433 2 : py /= rg->sfactor;
434 : }
435 132117 : PetscUseTypeMethod(rg,checkinside,px,py,inside+i);
436 132117 : if (PetscUnlikely(rg->complement)) inside[i] = -inside[i];
437 : }
438 66178 : PetscFunctionReturn(PETSC_SUCCESS);
439 : }
440 :
441 : /*@
442 : RGIsAxisymmetric - Determines if the region is symmetric with respect
443 : to the real or imaginary axis.
444 :
445 : Not Collective
446 :
447 : Input Parameters:
448 : + rg - the region context
449 : - vertical - true if symmetry must be checked against the vertical axis
450 :
451 : Output Parameter:
452 : . symm - true if the region is axisymmetric
453 :
454 : Note:
455 : If the vertical argument is true, symmetry is checked with respect to
456 : the vertical axis, otherwise with respect to the horizontal axis.
457 :
458 : Level: intermediate
459 :
460 : .seealso: RGCanUseConjugates()
461 : @*/
462 9 : PetscErrorCode RGIsAxisymmetric(RG rg,PetscBool vertical,PetscBool *symm)
463 : {
464 9 : PetscFunctionBegin;
465 9 : PetscValidHeaderSpecific(rg,RG_CLASSID,1);
466 9 : PetscValidType(rg,1);
467 9 : PetscAssertPointer(symm,3);
468 9 : *symm = PETSC_FALSE;
469 9 : PetscTryTypeMethod(rg,isaxisymmetric,vertical,symm);
470 9 : PetscFunctionReturn(PETSC_SUCCESS);
471 : }
472 :
473 : /*@
474 : RGCanUseConjugates - Used in contour integral methods to determine whether
475 : half of integration points can be avoided (use their conjugates).
476 :
477 : Not Collective
478 :
479 : Input Parameters:
480 : + rg - the region context
481 : - realmats - true if the problem matrices are real
482 :
483 : Output Parameter:
484 : . useconj - whether it is possible to use conjugates
485 :
486 : Notes:
487 : If some integration points are the conjugates of other points, then the
488 : associated computational cost can be saved. This depends on the problem
489 : matrices being real and also the region being symmetric with respect to
490 : the horizontal axis. The result is false if using real arithmetic or
491 : in the case of a flat region (height equal to zero).
492 :
493 : Level: developer
494 :
495 : .seealso: RGIsAxisymmetric()
496 : @*/
497 26 : PetscErrorCode RGCanUseConjugates(RG rg,PetscBool realmats,PetscBool *useconj)
498 : {
499 : #if defined(PETSC_USE_COMPLEX)
500 : PetscReal c,d;
501 : PetscBool isaxisymm;
502 : #endif
503 :
504 26 : PetscFunctionBegin;
505 26 : PetscValidHeaderSpecific(rg,RG_CLASSID,1);
506 26 : PetscValidType(rg,1);
507 26 : PetscAssertPointer(useconj,3);
508 26 : *useconj = PETSC_FALSE;
509 : #if defined(PETSC_USE_COMPLEX)
510 : if (realmats) {
511 : PetscCall(RGIsAxisymmetric(rg,PETSC_FALSE,&isaxisymm));
512 : if (isaxisymm) {
513 : PetscCall(RGComputeBoundingBox(rg,NULL,NULL,&c,&d));
514 : if (c!=d) *useconj = PETSC_TRUE;
515 : }
516 : }
517 : #endif
518 26 : PetscFunctionReturn(PETSC_SUCCESS);
519 : }
520 :
521 : /*@
522 : RGComputeContour - Computes the coordinates of several points lying on the
523 : contour of the region.
524 :
525 : Not Collective
526 :
527 : Input Parameters:
528 : + rg - the region context
529 : - n - number of points to compute
530 :
531 : Output Parameters:
532 : + cr - location to store real parts
533 : - ci - location to store imaginary parts
534 :
535 : Notes:
536 : In real scalars, either cr or ci can be NULL (but not both). In complex
537 : scalars, the coordinates are stored in cr, which cannot be NULL (ci is
538 : not referenced).
539 :
540 : Level: intermediate
541 :
542 : .seealso: RGComputeBoundingBox()
543 : @*/
544 64 : PetscErrorCode RGComputeContour(RG rg,PetscInt n,PetscScalar cr[],PetscScalar ci[])
545 : {
546 64 : PetscInt i;
547 :
548 64 : PetscFunctionBegin;
549 64 : PetscValidHeaderSpecific(rg,RG_CLASSID,1);
550 64 : PetscValidType(rg,1);
551 : #if defined(PETSC_USE_COMPLEX)
552 : PetscAssertPointer(cr,3);
553 : #else
554 64 : PetscCheck(cr || ci,PetscObjectComm((PetscObject)rg),PETSC_ERR_SUP,"cr and ci cannot be NULL simultaneously");
555 : #endif
556 64 : PetscCheck(!rg->complement,PetscObjectComm((PetscObject)rg),PETSC_ERR_SUP,"Cannot compute contour of region with complement flag set");
557 64 : PetscUseTypeMethod(rg,computecontour,n,cr,ci);
558 152159 : for (i=0;i<n;i++) {
559 152095 : if (cr) cr[i] *= rg->sfactor;
560 152095 : if (ci) ci[i] *= rg->sfactor;
561 : }
562 64 : PetscFunctionReturn(PETSC_SUCCESS);
563 : }
564 :
565 : /*@
566 : RGComputeBoundingBox - Determines the endpoints of a rectangle in the complex plane that
567 : contains the region.
568 :
569 : Not Collective
570 :
571 : Input Parameter:
572 : . rg - the region context
573 :
574 : Output Parameters:
575 : + a - left endpoint of the bounding box in the real axis
576 : . b - right endpoint of the bounding box in the real axis
577 : . c - bottom endpoint of the bounding box in the imaginary axis
578 : - d - top endpoint of the bounding box in the imaginary axis
579 :
580 : Notes:
581 : The bounding box is defined as [a,b]x[c,d]. In regions that are not bounded (e.g. an
582 : open interval) or with the complement flag set, it makes no sense to compute a bounding
583 : box, so the return values are infinite.
584 :
585 : Level: intermediate
586 :
587 : .seealso: RGComputeContour()
588 : @*/
589 13 : PetscErrorCode RGComputeBoundingBox(RG rg,PetscReal *a,PetscReal *b,PetscReal *c,PetscReal *d)
590 : {
591 13 : PetscFunctionBegin;
592 13 : PetscValidHeaderSpecific(rg,RG_CLASSID,1);
593 13 : PetscValidType(rg,1);
594 :
595 13 : if (rg->complement) { /* cannot compute bounding box */
596 0 : if (a) *a = -PETSC_MAX_REAL;
597 0 : if (b) *b = PETSC_MAX_REAL;
598 0 : if (c) *c = -PETSC_MAX_REAL;
599 0 : if (d) *d = PETSC_MAX_REAL;
600 : } else {
601 13 : PetscUseTypeMethod(rg,computebbox,a,b,c,d);
602 13 : if (a && *a!=-PETSC_MAX_REAL) *a *= rg->sfactor;
603 13 : if (b && *b!= PETSC_MAX_REAL) *b *= rg->sfactor;
604 13 : if (c && *c!=-PETSC_MAX_REAL) *c *= rg->sfactor;
605 13 : if (d && *d!= PETSC_MAX_REAL) *d *= rg->sfactor;
606 : }
607 13 : PetscFunctionReturn(PETSC_SUCCESS);
608 : }
609 :
610 : /*@
611 : RGComputeQuadrature - Computes the values of the parameters used in a
612 : quadrature rule for a contour integral around the boundary of the region.
613 :
614 : Not Collective
615 :
616 : Input Parameters:
617 : + rg - the region context
618 : . quad - the type of quadrature
619 : - n - number of quadrature points to compute
620 :
621 : Output Parameters:
622 : + z - quadrature points
623 : . zn - normalized quadrature points
624 : - w - quadrature weights
625 :
626 : Notes:
627 : In complex scalars, the values returned in z are often the same as those
628 : computed by RGComputeContour(), but this is not the case in real scalars
629 : where all output arguments are real.
630 :
631 : The computed values change for different quadrature rules (trapezoidal
632 : or Chebyshev).
633 :
634 : Level: intermediate
635 :
636 : .seealso: RGComputeContour()
637 : @*/
638 26 : PetscErrorCode RGComputeQuadrature(RG rg,RGQuadRule quad,PetscInt n,PetscScalar z[],PetscScalar zn[],PetscScalar w[])
639 : {
640 26 : PetscFunctionBegin;
641 26 : PetscValidHeaderSpecific(rg,RG_CLASSID,1);
642 26 : PetscValidType(rg,1);
643 26 : PetscAssertPointer(z,4);
644 26 : PetscAssertPointer(zn,5);
645 26 : PetscAssertPointer(w,6);
646 :
647 26 : PetscCall(RGComputeContour(rg,n,z,NULL));
648 26 : PetscUseTypeMethod(rg,computequadrature,quad,n,z,zn,w);
649 26 : PetscFunctionReturn(PETSC_SUCCESS);
650 : }
651 :
652 : /*@
653 : RGSetComplement - Sets a flag to indicate that the region is the complement
654 : of the specified one.
655 :
656 : Logically Collective
657 :
658 : Input Parameters:
659 : + rg - the region context
660 : - flg - the boolean flag
661 :
662 : Options Database Key:
663 : . -rg_complement <bool> - Activate/deactivate the complementation of the region
664 :
665 : Level: intermediate
666 :
667 : .seealso: RGGetComplement()
668 : @*/
669 1 : PetscErrorCode RGSetComplement(RG rg,PetscBool flg)
670 : {
671 1 : PetscFunctionBegin;
672 1 : PetscValidHeaderSpecific(rg,RG_CLASSID,1);
673 3 : PetscValidLogicalCollectiveBool(rg,flg,2);
674 1 : rg->complement = flg;
675 1 : PetscFunctionReturn(PETSC_SUCCESS);
676 : }
677 :
678 : /*@
679 : RGGetComplement - Gets a flag that indicates whether the region
680 : is complemented or not.
681 :
682 : Not Collective
683 :
684 : Input Parameter:
685 : . rg - the region context
686 :
687 : Output Parameter:
688 : . flg - the flag
689 :
690 : Level: intermediate
691 :
692 : .seealso: RGSetComplement()
693 : @*/
694 26 : PetscErrorCode RGGetComplement(RG rg,PetscBool *flg)
695 : {
696 26 : PetscFunctionBegin;
697 26 : PetscValidHeaderSpecific(rg,RG_CLASSID,1);
698 26 : PetscAssertPointer(flg,2);
699 26 : *flg = rg->complement;
700 26 : PetscFunctionReturn(PETSC_SUCCESS);
701 : }
702 :
703 : /*@
704 : RGSetScale - Sets the scaling factor to be used when checking that a
705 : point is inside the region and when computing the contour.
706 :
707 : Logically Collective
708 :
709 : Input Parameters:
710 : + rg - the region context
711 : - sfactor - the scaling factor
712 :
713 : Options Database Key:
714 : . -rg_scale <real> - Sets the scaling factor
715 :
716 : Level: advanced
717 :
718 : .seealso: RGGetScale(), RGCheckInside()
719 : @*/
720 1 : PetscErrorCode RGSetScale(RG rg,PetscReal sfactor)
721 : {
722 1 : PetscFunctionBegin;
723 1 : PetscValidHeaderSpecific(rg,RG_CLASSID,1);
724 3 : PetscValidLogicalCollectiveReal(rg,sfactor,2);
725 1 : if (sfactor == (PetscReal)PETSC_DEFAULT || sfactor == (PetscReal)PETSC_DECIDE) sfactor = 1.0;
726 1 : PetscCheck(sfactor>0.0,PetscObjectComm((PetscObject)rg),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of scaling factor. Must be > 0");
727 1 : rg->sfactor = sfactor;
728 1 : PetscFunctionReturn(PETSC_SUCCESS);
729 : }
730 :
731 : /*@
732 : RGGetScale - Gets the scaling factor.
733 :
734 : Not Collective
735 :
736 : Input Parameter:
737 : . rg - the region context
738 :
739 : Output Parameter:
740 : . sfactor - the scaling factor
741 :
742 : Level: advanced
743 :
744 : .seealso: RGSetScale()
745 : @*/
746 52 : PetscErrorCode RGGetScale(RG rg,PetscReal *sfactor)
747 : {
748 52 : PetscFunctionBegin;
749 52 : PetscValidHeaderSpecific(rg,RG_CLASSID,1);
750 52 : PetscAssertPointer(sfactor,2);
751 52 : *sfactor = rg->sfactor;
752 52 : PetscFunctionReturn(PETSC_SUCCESS);
753 : }
754 :
755 : /*@
756 : RGPushScale - Sets an additional scaling factor, that will multiply the
757 : user-defined scaling factor.
758 :
759 : Logically Collective
760 :
761 : Input Parameters:
762 : + rg - the region context
763 : - sfactor - the scaling factor
764 :
765 : Notes:
766 : The current implementation does not allow pushing several scaling factors.
767 :
768 : This is intended for internal use, for instance in polynomial eigensolvers
769 : that use parameter scaling.
770 :
771 : Level: developer
772 :
773 : .seealso: RGPopScale(), RGSetScale()
774 : @*/
775 79 : PetscErrorCode RGPushScale(RG rg,PetscReal sfactor)
776 : {
777 79 : PetscFunctionBegin;
778 79 : PetscValidHeaderSpecific(rg,RG_CLASSID,1);
779 237 : PetscValidLogicalCollectiveReal(rg,sfactor,2);
780 79 : PetscCheck(sfactor>0.0,PetscObjectComm((PetscObject)rg),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of scaling factor. Must be > 0");
781 79 : PetscCheck(!rg->osfactor,PetscObjectComm((PetscObject)rg),PETSC_ERR_SUP,"Current implementation does not allow pushing several scaling factors");
782 79 : rg->osfactor = rg->sfactor;
783 79 : rg->sfactor *= sfactor;
784 79 : PetscFunctionReturn(PETSC_SUCCESS);
785 : }
786 :
787 : /*@
788 : RGPopScale - Pops the scaling factor set with RGPushScale().
789 :
790 : Logically Collective
791 :
792 : Input Parameter:
793 : . rg - the region context
794 :
795 : Level: developer
796 :
797 : .seealso: RGPushScale()
798 : @*/
799 79 : PetscErrorCode RGPopScale(RG rg)
800 : {
801 79 : PetscFunctionBegin;
802 79 : PetscValidHeaderSpecific(rg,RG_CLASSID,1);
803 79 : PetscCheck(rg->osfactor,PetscObjectComm((PetscObject)rg),PETSC_ERR_ORDER,"Must call RGPushScale first");
804 79 : rg->sfactor = rg->osfactor;
805 79 : rg->osfactor = 0.0;
806 79 : PetscFunctionReturn(PETSC_SUCCESS);
807 : }
808 :
809 : /*@
810 : RGDestroy - Destroys RG context that was created with RGCreate().
811 :
812 : Collective
813 :
814 : Input Parameter:
815 : . rg - the region context
816 :
817 : Level: beginner
818 :
819 : .seealso: RGCreate()
820 : @*/
821 954 : PetscErrorCode RGDestroy(RG *rg)
822 : {
823 954 : PetscFunctionBegin;
824 954 : if (!*rg) PetscFunctionReturn(PETSC_SUCCESS);
825 930 : PetscValidHeaderSpecific(*rg,RG_CLASSID,1);
826 930 : if (--((PetscObject)*rg)->refct > 0) { *rg = NULL; PetscFunctionReturn(PETSC_SUCCESS); }
827 920 : PetscTryTypeMethod(*rg,destroy);
828 920 : PetscCall(PetscHeaderDestroy(rg));
829 920 : PetscFunctionReturn(PETSC_SUCCESS);
830 : }
831 :
832 : /*@C
833 : RGRegister - Adds a region to the RG package.
834 :
835 : Not Collective
836 :
837 : Input Parameters:
838 : + name - name of a new user-defined RG
839 : - function - routine to create context
840 :
841 : Notes:
842 : RGRegister() may be called multiple times to add several user-defined regions.
843 :
844 : Level: advanced
845 :
846 : .seealso: RGRegisterAll()
847 : @*/
848 3144 : PetscErrorCode RGRegister(const char *name,PetscErrorCode (*function)(RG))
849 : {
850 3144 : PetscFunctionBegin;
851 3144 : PetscCall(RGInitializePackage());
852 3144 : PetscCall(PetscFunctionListAdd(&RGList,name,function));
853 3144 : PetscFunctionReturn(PETSC_SUCCESS);
854 : }
|
