Actual source code: dmproject.c
petsc-3.7.7 2017-09-25
2: #include <petsc/private/petscimpl.h>
3: #include <petscdm.h> /*I "petscdm.h" I*/
4: #include <petscdmplex.h> /*I "petscdmplex.h" I*/
5: #include <petscksp.h> /*I "petscksp.h" I*/
7: typedef struct _projectConstraintsCtx
8: {
9: DM dm;
10: Vec mask;
11: }
12: projectConstraintsCtx;
16: PetscErrorCode MatMult_GlobalToLocalNormal(Mat CtC, Vec x, Vec y)
17: {
18: DM dm;
19: Vec local, mask;
20: projectConstraintsCtx *ctx;
21: PetscErrorCode ierr;
24: MatShellGetContext(CtC,&ctx);
25: dm = ctx->dm;
26: mask = ctx->mask;
27: DMGetLocalVector(dm,&local);
28: DMGlobalToLocalBegin(dm,x,INSERT_VALUES,local);
29: DMGlobalToLocalEnd(dm,x,INSERT_VALUES,local);
30: if (mask) {VecPointwiseMult(local,mask,local);}
31: VecSet(y,0.);
32: DMLocalToGlobalBegin(dm,local,ADD_VALUES,y);
33: DMLocalToGlobalEnd(dm,local,ADD_VALUES,y);
34: DMRestoreLocalVector(dm,&local);
35: return(0);
36: }
40: static PetscErrorCode DMGlobalToLocalSolve_project1 (PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar u[], void *ctx)
41: {
42: PetscInt f;
45: for (f = 0; f < Nf; f++) {
46: u[f] = 1.;
47: }
48: return(0);
49: }
53: /*@
54: DMGlobalToLocalSolve - Solve for the global vector that is mapped to a given local vector by DMGlobalToLocalBegin()/DMGlobalToLocalEnd() with mode
55: = INSERT_VALUES. It is assumed that the sum of all the local vector sizes is greater than or equal to the global vector size, so the solution is
56: a least-squares solution. It is also assumed that DMLocalToGlobalBegin()/DMLocalToGlobalEnd() with mode = ADD_VALUES is the adjoint of the
57: global-to-local map, so that the least-squares solution may be found by the normal equations.
59: collective
61: Input Parameters:
62: + dm - The DM object
63: . x - The local vector
64: - y - The global vector: the input value of globalVec is used as an initial guess
66: Output Parameters:
67: . y - The least-squares solution
69: Level: advanced
71: Note: If the DM is of type DMPLEX, then y is the solution of L' * D * L * y = L' * D * x, where D is a diagonal mask that is 1 for every point in
72: the union of the closures of the local cells and 0 otherwise. This difference is only relevant if there are anchor points that are not in the
73: closure of any local cell (see DMPlexGetAnchors()/DMPlexSetAnchors()).
75: .seealso: DMGlobalToLocalBegin(), DMGlobalToLocalEnd(), DMLocalToGlobalBegin(), DMLocalToGlobalEnd(), DMPlexGetAnchors(), DMPlexSetAnchors()
76: @*/
77: PetscErrorCode DMGlobalToLocalSolve(DM dm, Vec x, Vec y)
78: {
79: Mat CtC;
80: PetscInt n, N, cStart, cEnd, cEndInterior, c;
81: PetscBool isPlex;
82: KSP ksp;
83: PC pc;
84: Vec global, mask=NULL;
85: projectConstraintsCtx ctx;
89: PetscObjectTypeCompare((PetscObject)dm,DMPLEX,&isPlex);
90: if (isPlex) {
91: /* mark points in the closure */
92: DMCreateLocalVector(dm,&mask);
93: VecSet(mask,0.0);
94: DMPlexGetHeightStratum(dm,0,&cStart,&cEnd);
95: DMPlexGetHybridBounds(dm,&cEndInterior,NULL,NULL,NULL);
96: cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
97: if (cEnd > cStart) {
98: PetscScalar *ones;
99: PetscInt numValues, i;
101: DMPlexVecGetClosure(dm,NULL,mask,cStart,&numValues,NULL);
102: PetscMalloc1(numValues,&ones);
103: for (i = 0; i < numValues; i++) {
104: ones[i] = 1.;
105: }
106: for (c = cStart; c < cEnd; c++) {
107: DMPlexVecSetClosure(dm,NULL,mask,c,ones,INSERT_VALUES);
108: }
109: PetscFree(ones);
110: }
111: }
112: else {
113: PetscBool hasMask;
115: DMHasNamedLocalVector(dm, "_DMGlobalToLocalSolve_mask", &hasMask);
116: if (!hasMask) {
117: PetscErrorCode (**func) (PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx);
118: void **ctx;
119: PetscInt Nf, f;
121: DMGetNumFields(dm, &Nf);
122: PetscMalloc2(Nf, &func, Nf, &ctx);
123: for (f = 0; f < Nf; ++f) {
124: func[f] = DMGlobalToLocalSolve_project1;
125: ctx[f] = NULL;
126: }
127: DMGetNamedLocalVector(dm, "_DMGlobalToLocalSolve_mask", &mask);
128: DMProjectFunctionLocal(dm,0.0,func,ctx,INSERT_ALL_VALUES,mask);
129: DMRestoreNamedLocalVector(dm, "_DMGlobalToLocalSolve_mask", &mask);
130: PetscFree2(func, ctx);
131: }
132: DMGetNamedLocalVector(dm, "_DMGlobalToLocalSolve_mask", &mask);
133: }
134: ctx.dm = dm;
135: ctx.mask = mask;
136: VecGetSize(y,&N);
137: VecGetLocalSize(y,&n);
138: MatCreate(PetscObjectComm((PetscObject)dm),&CtC);
139: MatSetSizes(CtC,n,n,N,N);
140: MatSetType(CtC,MATSHELL);
141: MatSetUp(CtC);
142: MatShellSetContext(CtC,&ctx);
143: MatShellSetOperation(CtC,MATOP_MULT,(void(*)(void))MatMult_GlobalToLocalNormal);
144: KSPCreate(PetscObjectComm((PetscObject)dm),&ksp);
145: KSPSetOperators(ksp,CtC,CtC);
146: KSPSetType(ksp,KSPCG);
147: KSPGetPC(ksp,&pc);
148: PCSetType(pc,PCNONE);
149: KSPSetInitialGuessNonzero(ksp,PETSC_TRUE);
150: KSPSetUp(ksp);
151: DMGetGlobalVector(dm,&global);
152: VecSet(global,0.);
153: if (mask) {VecPointwiseMult(x,mask,x);}
154: DMLocalToGlobalBegin(dm,x,ADD_VALUES,global);
155: DMLocalToGlobalEnd(dm,x,ADD_VALUES,global);
156: KSPSolve(ksp,global,y);
157: DMRestoreGlobalVector(dm,&global);
158: /* clean up */
159: KSPDestroy(&ksp);
160: MatDestroy(&CtC);
161: if (isPlex) {
162: VecDestroy(&mask);
163: }
164: else {
165: DMRestoreNamedLocalVector(dm, "_DMGlobalToLocalSolve_mask", &mask);
166: }
168: return(0);
169: }
173: /*@C
174: DMProjectField - This projects the given function of the fields into the function space provided.
176: Input Parameters:
177: + dm - The DM
178: . U - The input field vector
179: . funcs - The functions to evaluate, one per field
180: - mode - The insertion mode for values
182: Output Parameter:
183: . X - The output vector
185: Level: developer
187: .seealso: DMProjectFunction(), DMComputeL2Diff()
188: @*/
189: PetscErrorCode DMProjectField(DM dm, Vec U,
190: void (**funcs)(PetscInt, PetscInt, PetscInt,
191: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
192: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
193: PetscReal, const PetscReal[], PetscScalar[]),
194: InsertMode mode, Vec X)
195: {
196: Vec localX, localU;
201: DMGetLocalVector(dm, &localX);
202: DMGetLocalVector(dm, &localU);
203: DMGlobalToLocalBegin(dm, U, INSERT_VALUES, localU);
204: DMGlobalToLocalEnd(dm, U, INSERT_VALUES, localU);
205: DMProjectFieldLocal(dm, localU, funcs, mode, localX);
206: DMLocalToGlobalBegin(dm, localX, mode, X);
207: DMLocalToGlobalEnd(dm, localX, mode, X);
208: if (mode == INSERT_VALUES || mode == INSERT_ALL_VALUES || mode == INSERT_BC_VALUES) {
209: Mat cMat;
211: DMGetDefaultConstraints(dm, NULL, &cMat);
212: if (cMat) {
213: DMGlobalToLocalSolve(dm, localX, X);
214: }
215: }
216: DMRestoreLocalVector(dm, &localX);
217: DMRestoreLocalVector(dm, &localU);
218: return(0);
219: }