Actual source code: dmproject.c
petsc-3.8.4 2018-03-24
2: #include <petsc/private/petscimpl.h>
3: #include <petscdm.h>
4: #include <petscdmplex.h>
5: #include <petscksp.h>
7: typedef struct _projectConstraintsCtx
8: {
9: DM dm;
10: Vec mask;
11: }
12: projectConstraintsCtx;
14: PetscErrorCode MatMult_GlobalToLocalNormal(Mat CtC, Vec x, Vec y)
15: {
16: DM dm;
17: Vec local, mask;
18: projectConstraintsCtx *ctx;
19: PetscErrorCode ierr;
22: MatShellGetContext(CtC,&ctx);
23: dm = ctx->dm;
24: mask = ctx->mask;
25: DMGetLocalVector(dm,&local);
26: DMGlobalToLocalBegin(dm,x,INSERT_VALUES,local);
27: DMGlobalToLocalEnd(dm,x,INSERT_VALUES,local);
28: if (mask) {VecPointwiseMult(local,mask,local);}
29: VecSet(y,0.);
30: DMLocalToGlobalBegin(dm,local,ADD_VALUES,y);
31: DMLocalToGlobalEnd(dm,local,ADD_VALUES,y);
32: DMRestoreLocalVector(dm,&local);
33: return(0);
34: }
36: static PetscErrorCode DMGlobalToLocalSolve_project1 (PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar u[], void *ctx)
37: {
38: PetscInt f;
41: for (f = 0; f < Nf; f++) {
42: u[f] = 1.;
43: }
44: return(0);
45: }
47: /*@
48: DMGlobalToLocalSolve - Solve for the global vector that is mapped to a given local vector by DMGlobalToLocalBegin()/DMGlobalToLocalEnd() with mode
49: = 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
50: a least-squares solution. It is also assumed that DMLocalToGlobalBegin()/DMLocalToGlobalEnd() with mode = ADD_VALUES is the adjoint of the
51: global-to-local map, so that the least-squares solution may be found by the normal equations.
53: collective
55: Input Parameters:
56: + dm - The DM object
57: . x - The local vector
58: - y - The global vector: the input value of globalVec is used as an initial guess
60: Output Parameters:
61: . y - The least-squares solution
63: Level: advanced
65: 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
66: 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
67: closure of any local cell (see DMPlexGetAnchors()/DMPlexSetAnchors()).
69: .seealso: DMGlobalToLocalBegin(), DMGlobalToLocalEnd(), DMLocalToGlobalBegin(), DMLocalToGlobalEnd(), DMPlexGetAnchors(), DMPlexSetAnchors()
70: @*/
71: PetscErrorCode DMGlobalToLocalSolve(DM dm, Vec x, Vec y)
72: {
73: Mat CtC;
74: PetscInt n, N, cStart, cEnd, cEndInterior, c;
75: PetscBool isPlex;
76: KSP ksp;
77: PC pc;
78: Vec global, mask=NULL;
79: projectConstraintsCtx ctx;
83: PetscObjectTypeCompare((PetscObject)dm,DMPLEX,&isPlex);
84: if (isPlex) {
85: /* mark points in the closure */
86: DMCreateLocalVector(dm,&mask);
87: VecSet(mask,0.0);
88: DMPlexGetHeightStratum(dm,0,&cStart,&cEnd);
89: DMPlexGetHybridBounds(dm,&cEndInterior,NULL,NULL,NULL);
90: cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
91: if (cEnd > cStart) {
92: PetscScalar *ones;
93: PetscInt numValues, i;
95: DMPlexVecGetClosure(dm,NULL,mask,cStart,&numValues,NULL);
96: PetscMalloc1(numValues,&ones);
97: for (i = 0; i < numValues; i++) {
98: ones[i] = 1.;
99: }
100: for (c = cStart; c < cEnd; c++) {
101: DMPlexVecSetClosure(dm,NULL,mask,c,ones,INSERT_VALUES);
102: }
103: PetscFree(ones);
104: }
105: }
106: else {
107: PetscBool hasMask;
109: DMHasNamedLocalVector(dm, "_DMGlobalToLocalSolve_mask", &hasMask);
110: if (!hasMask) {
111: PetscErrorCode (**func) (PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx);
112: void **ctx;
113: PetscInt Nf, f;
115: DMGetNumFields(dm, &Nf);
116: PetscMalloc2(Nf, &func, Nf, &ctx);
117: for (f = 0; f < Nf; ++f) {
118: func[f] = DMGlobalToLocalSolve_project1;
119: ctx[f] = NULL;
120: }
121: DMGetNamedLocalVector(dm, "_DMGlobalToLocalSolve_mask", &mask);
122: DMProjectFunctionLocal(dm,0.0,func,ctx,INSERT_ALL_VALUES,mask);
123: DMRestoreNamedLocalVector(dm, "_DMGlobalToLocalSolve_mask", &mask);
124: PetscFree2(func, ctx);
125: }
126: DMGetNamedLocalVector(dm, "_DMGlobalToLocalSolve_mask", &mask);
127: }
128: ctx.dm = dm;
129: ctx.mask = mask;
130: VecGetSize(y,&N);
131: VecGetLocalSize(y,&n);
132: MatCreate(PetscObjectComm((PetscObject)dm),&CtC);
133: MatSetSizes(CtC,n,n,N,N);
134: MatSetType(CtC,MATSHELL);
135: MatSetUp(CtC);
136: MatShellSetContext(CtC,&ctx);
137: MatShellSetOperation(CtC,MATOP_MULT,(void(*)(void))MatMult_GlobalToLocalNormal);
138: KSPCreate(PetscObjectComm((PetscObject)dm),&ksp);
139: KSPSetOperators(ksp,CtC,CtC);
140: KSPSetType(ksp,KSPCG);
141: KSPGetPC(ksp,&pc);
142: PCSetType(pc,PCNONE);
143: KSPSetInitialGuessNonzero(ksp,PETSC_TRUE);
144: KSPSetUp(ksp);
145: DMGetGlobalVector(dm,&global);
146: VecSet(global,0.);
147: if (mask) {VecPointwiseMult(x,mask,x);}
148: DMLocalToGlobalBegin(dm,x,ADD_VALUES,global);
149: DMLocalToGlobalEnd(dm,x,ADD_VALUES,global);
150: KSPSolve(ksp,global,y);
151: DMRestoreGlobalVector(dm,&global);
152: /* clean up */
153: KSPDestroy(&ksp);
154: MatDestroy(&CtC);
155: if (isPlex) {
156: VecDestroy(&mask);
157: }
158: else {
159: DMRestoreNamedLocalVector(dm, "_DMGlobalToLocalSolve_mask", &mask);
160: }
162: return(0);
163: }
165: /*@C
166: DMProjectField - This projects the given function of the fields into the function space provided.
168: Input Parameters:
169: + dm - The DM
170: . time - The time
171: . U - The input field vector
172: . funcs - The functions to evaluate, one per field
173: - mode - The insertion mode for values
175: Output Parameter:
176: . X - The output vector
178: Level: developer
180: .seealso: DMProjectFunction(), DMComputeL2Diff()
181: @*/
182: PetscErrorCode DMProjectField(DM dm, PetscReal time, Vec U,
183: void (**funcs)(PetscInt, PetscInt, PetscInt,
184: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
185: const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
186: PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]),
187: InsertMode mode, Vec X)
188: {
189: Vec localX, localU;
194: DMGetLocalVector(dm, &localX);
195: /* We currently check whether locU == locX to see if we need to apply BC */
196: if (U != X) {DMGetLocalVector(dm, &localU);}
197: else {localU = localX;}
198: DMGlobalToLocalBegin(dm, U, INSERT_VALUES, localU);
199: DMGlobalToLocalEnd(dm, U, INSERT_VALUES, localU);
200: DMProjectFieldLocal(dm, time, localU, funcs, mode, localX);
201: DMLocalToGlobalBegin(dm, localX, mode, X);
202: DMLocalToGlobalEnd(dm, localX, mode, X);
203: if (mode == INSERT_VALUES || mode == INSERT_ALL_VALUES || mode == INSERT_BC_VALUES) {
204: Mat cMat;
206: DMGetDefaultConstraints(dm, NULL, &cMat);
207: if (cMat) {
208: DMGlobalToLocalSolve(dm, localX, X);
209: }
210: }
211: DMRestoreLocalVector(dm, &localX);
212: if (U != X) {DMRestoreLocalVector(dm, &localU);}
213: return(0);
214: }