1 #include <petsc/private/dmfieldimpl.h> /*I "petscdmfield.h" I*/
2 #include <petsc/private/dmimpl.h> /*I "petscdm.h" I*/
3 #include <petscdmda.h>
4
5 typedef struct _n_DMField_DA {
6 PetscScalar *cornerVals;
7 PetscScalar *cornerCoeffs;
8 PetscScalar *work;
9 PetscReal coordRange[3][2];
10 } DMField_DA;
11
DMFieldDestroy_DA(DMField field)12 static PetscErrorCode DMFieldDestroy_DA(DMField field)
13 {
14 DMField_DA *dafield;
15
16 PetscFunctionBegin;
17 dafield = (DMField_DA *)field->data;
18 PetscCall(PetscFree3(dafield->cornerVals, dafield->cornerCoeffs, dafield->work));
19 PetscCall(PetscFree(dafield));
20 PetscFunctionReturn(PETSC_SUCCESS);
21 }
22
DMFieldView_DA(DMField field,PetscViewer viewer)23 static PetscErrorCode DMFieldView_DA(DMField field, PetscViewer viewer)
24 {
25 DMField_DA *dafield = (DMField_DA *)field->data;
26 PetscBool isascii;
27
28 PetscFunctionBegin;
29 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
30 if (isascii) {
31 PetscInt i, dim;
32 PetscInt nc;
33 DM dm = field->dm;
34
35 PetscCall(PetscViewerASCIIPrintf(viewer, "Field corner values:\n"));
36 PetscCall(PetscViewerASCIIPushTab(viewer));
37 PetscCall(DMGetDimension(dm, &dim));
38 nc = field->numComponents;
39 for (i = 0; i < (1 << dim); i++) {
40 PetscInt j;
41
42 for (j = 0; j < nc; j++) {
43 PetscScalar val = dafield->cornerVals[nc * i + j];
44
45 #if !defined(PETSC_USE_COMPLEX)
46 PetscCall(PetscViewerASCIIPrintf(viewer, "%g ", (double)val));
47 #else
48 PetscCall(PetscViewerASCIIPrintf(viewer, "%g+i%g ", (double)PetscRealPart(val), (double)PetscImaginaryPart(val)));
49 #endif
50 }
51 PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
52 }
53 PetscCall(PetscViewerASCIIPopTab(viewer));
54 }
55 PetscFunctionReturn(PETSC_SUCCESS);
56 }
57
58 #define MEdot(y, A, x, m, c, cast) \
59 do { \
60 PetscInt _k, _l; \
61 for (_k = 0; _k < (c); _k++) (y)[_k] = 0.; \
62 for (_l = 0; _l < (m); _l++) { \
63 for (_k = 0; _k < (c); _k++) (y)[_k] += cast((A)[(c) * _l + _k]) * (x)[_l]; \
64 } \
65 } while (0)
66
67 #define MEHess(out, cf, etaB, etaD, dim, nc, cast) \
68 do { \
69 PetscInt _m, _j, _k; \
70 for (_m = 0; _m < (nc) * (dim) * (dim); _m++) (out)[_m] = 0.; \
71 for (_j = 0; _j < (dim); _j++) { \
72 for (_k = _j + 1; _k < (dim); _k++) { \
73 PetscInt _ind = (1 << _j) + (1 << _k); \
74 for (_m = 0; _m < (nc); _m++) { \
75 PetscScalar c = (cf)[_m] * (etaB)[_ind] * (etaD)[_ind]; \
76 (out)[(_m * (dim) + _k) * (dim) + _j] += cast(c); \
77 (out)[(_m * (dim) + _j) * (dim) + _k] += cast(c); \
78 } \
79 } \
80 } \
81 } while (0)
82
MultilinearEvaluate(PetscInt dim,PetscReal (* coordRange)[2],PetscInt nc,PetscScalar * cf,PetscScalar * cfWork,PetscInt nPoints,const PetscScalar * points,PetscDataType datatype,void * B,void * D,void * H)83 static void MultilinearEvaluate(PetscInt dim, PetscReal (*coordRange)[2], PetscInt nc, PetscScalar *cf, PetscScalar *cfWork, PetscInt nPoints, const PetscScalar *points, PetscDataType datatype, void *B, void *D, void *H)
84 {
85 PetscInt i, j, k, l, m;
86 PetscInt whol = 1 << dim;
87 PetscInt half = whol >> 1;
88
89 PetscFunctionBeginHot;
90 if (!B && !D && !H) PetscFunctionReturnVoid();
91 for (i = 0; i < nPoints; i++) {
92 const PetscScalar *point = &points[dim * i];
93 PetscReal deta[3] = {0.};
94 PetscReal etaB[8] = {1., 1., 1., 1., 1., 1., 1., 1.};
95 PetscReal etaD[8] = {1., 1., 1., 1., 1., 1., 1., 1.};
96 PetscReal work[8];
97
98 for (j = 0; j < dim; j++) {
99 PetscReal e, d;
100
101 e = (PetscRealPart(point[j]) - coordRange[j][0]) / coordRange[j][1];
102 deta[j] = d = 1. / coordRange[j][1];
103 for (k = 0; k < whol; k++) work[k] = etaB[k];
104 for (k = 0; k < half; k++) {
105 etaB[k] = work[2 * k] * e;
106 etaB[k + half] = work[2 * k + 1];
107 }
108 if (H) {
109 for (k = 0; k < whol; k++) work[k] = etaD[k];
110 for (k = 0; k < half; k++) {
111 etaD[k + half] = work[2 * k];
112 etaD[k] = work[2 * k + 1] * d;
113 }
114 }
115 }
116 if (B) {
117 if (datatype == PETSC_SCALAR) {
118 PetscScalar *out = &((PetscScalar *)B)[nc * i];
119
120 MEdot(out, cf, etaB, (1 << dim), nc, (PetscScalar));
121 } else {
122 PetscReal *out = &((PetscReal *)B)[nc * i];
123
124 MEdot(out, cf, etaB, (1 << dim), nc, PetscRealPart);
125 }
126 }
127 if (D) {
128 if (datatype == PETSC_SCALAR) {
129 PetscScalar *out = &((PetscScalar *)D)[nc * dim * i];
130
131 for (m = 0; m < nc * dim; m++) out[m] = 0.;
132 } else {
133 PetscReal *out = &((PetscReal *)D)[nc * dim * i];
134
135 for (m = 0; m < nc * dim; m++) out[m] = 0.;
136 }
137 for (j = 0; j < dim; j++) {
138 PetscReal d = deta[j];
139
140 for (k = 0; k < whol * nc; k++) cfWork[k] = cf[k];
141 for (k = 0; k < whol; k++) work[k] = etaB[k];
142 for (k = 0; k < half; k++) {
143 PetscReal e;
144
145 etaB[k] = work[2 * k];
146 etaB[k + half] = e = work[2 * k + 1];
147
148 for (l = 0; l < nc; l++) {
149 cf[k * nc + l] = cfWork[2 * k * nc + l];
150 cf[(k + half) * nc + l] = cfWork[(2 * k + 1) * nc + l];
151 }
152 if (datatype == PETSC_SCALAR) {
153 PetscScalar *out = &((PetscScalar *)D)[nc * dim * i];
154
155 for (l = 0; l < nc; l++) out[l * dim + j] += d * e * cf[k * nc + l];
156 } else {
157 PetscReal *out = &((PetscReal *)D)[nc * dim * i];
158
159 for (l = 0; l < nc; l++) out[l * dim + j] += d * e * PetscRealPart(cf[k * nc + l]);
160 }
161 }
162 }
163 }
164 if (H) {
165 if (datatype == PETSC_SCALAR) {
166 PetscScalar *out = &((PetscScalar *)H)[nc * dim * dim * i];
167
168 MEHess(out, cf, etaB, etaD, dim, nc, (PetscScalar));
169 } else {
170 PetscReal *out = &((PetscReal *)H)[nc * dim * dim * i];
171
172 MEHess(out, cf, etaB, etaD, dim, nc, PetscRealPart);
173 }
174 }
175 }
176 PetscFunctionReturnVoid();
177 }
178
DMFieldEvaluate_DA(DMField field,Vec points,PetscDataType datatype,void * B,void * D,void * H)179 static PetscErrorCode DMFieldEvaluate_DA(DMField field, Vec points, PetscDataType datatype, void *B, void *D, void *H)
180 {
181 DM dm;
182 DMField_DA *dafield;
183 PetscInt dim;
184 PetscInt N, n, nc;
185 const PetscScalar *array;
186 PetscReal (*coordRange)[2];
187
188 PetscFunctionBegin;
189 dm = field->dm;
190 nc = field->numComponents;
191 dafield = (DMField_DA *)field->data;
192 PetscCall(DMGetDimension(dm, &dim));
193 PetscCall(VecGetLocalSize(points, &N));
194 PetscCheck(N % dim == 0, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Point vector size %" PetscInt_FMT " not divisible by coordinate dimension %" PetscInt_FMT, N, dim);
195 n = N / dim;
196 coordRange = &dafield->coordRange[0];
197 PetscCall(VecGetArrayRead(points, &array));
198 MultilinearEvaluate(dim, coordRange, nc, dafield->cornerCoeffs, dafield->work, n, array, datatype, B, D, H);
199 PetscCall(VecRestoreArrayRead(points, &array));
200 PetscFunctionReturn(PETSC_SUCCESS);
201 }
202
DMFieldEvaluateFE_DA(DMField field,IS cellIS,PetscQuadrature points,PetscDataType datatype,void * B,void * D,void * H)203 static PetscErrorCode DMFieldEvaluateFE_DA(DMField field, IS cellIS, PetscQuadrature points, PetscDataType datatype, void *B, void *D, void *H)
204 {
205 PetscInt c, i, j, k, dim, cellsPer[3] = {0}, first[3] = {0}, whol, half;
206 PetscReal stepPer[3] = {0.};
207 PetscReal cellCoordRange[3][2] = {
208 {0., 1.},
209 {0., 1.},
210 {0., 1.}
211 };
212 PetscScalar *cellCoeffs, *work;
213 DM dm;
214 DMDALocalInfo info;
215 PetscInt cStart, cEnd;
216 PetscInt nq, nc;
217 const PetscReal *q;
218 #if defined(PETSC_USE_COMPLEX)
219 PetscScalar *qs;
220 #else
221 const PetscScalar *qs;
222 #endif
223 DMField_DA *dafield;
224 PetscBool isStride;
225 const PetscInt *cells = NULL;
226 PetscInt sfirst = -1, stride = -1, nCells;
227
228 PetscFunctionBegin;
229 dafield = (DMField_DA *)field->data;
230 dm = field->dm;
231 nc = field->numComponents;
232 PetscCall(DMDAGetLocalInfo(dm, &info));
233 dim = info.dim;
234 work = dafield->work;
235 stepPer[0] = 1. / info.mx;
236 stepPer[1] = 1. / info.my;
237 stepPer[2] = 1. / info.mz;
238 first[0] = info.gxs;
239 first[1] = info.gys;
240 first[2] = info.gzs;
241 cellsPer[0] = info.gxm;
242 cellsPer[1] = info.gym;
243 cellsPer[2] = info.gzm;
244 /* TODO: probably take components into account */
245 PetscCall(PetscQuadratureGetData(points, NULL, NULL, &nq, &q, NULL));
246 #if defined(PETSC_USE_COMPLEX)
247 PetscCall(DMGetWorkArray(dm, nq * dim, MPIU_SCALAR, &qs));
248 for (i = 0; i < nq * dim; i++) qs[i] = q[i];
249 #else
250 qs = q;
251 #endif
252 PetscCall(DMDAGetHeightStratum(dm, 0, &cStart, &cEnd));
253 PetscCall(DMGetWorkArray(dm, (1 << dim) * nc, MPIU_SCALAR, &cellCoeffs));
254 whol = (1 << dim);
255 half = whol >> 1;
256 PetscCall(ISGetLocalSize(cellIS, &nCells));
257 PetscCall(PetscObjectTypeCompare((PetscObject)cellIS, ISSTRIDE, &isStride));
258 if (isStride) PetscCall(ISStrideGetInfo(cellIS, &sfirst, &stride));
259 else PetscCall(ISGetIndices(cellIS, &cells));
260 for (c = 0; c < nCells; c++) {
261 PetscInt cell = isStride ? (sfirst + c * stride) : cells[c];
262 PetscInt rem = cell;
263 PetscInt ijk[3] = {0};
264 void *cB, *cD, *cH;
265
266 if (datatype == PETSC_SCALAR) {
267 cB = PetscSafePointerPlusOffset((PetscScalar *)B, nc * nq * c);
268 cD = PetscSafePointerPlusOffset((PetscScalar *)D, nc * nq * dim * c);
269 cH = PetscSafePointerPlusOffset((PetscScalar *)H, nc * nq * dim * dim * c);
270 } else {
271 cB = PetscSafePointerPlusOffset((PetscReal *)B, nc * nq * c);
272 cD = PetscSafePointerPlusOffset((PetscReal *)D, nc * nq * dim * c);
273 cH = PetscSafePointerPlusOffset((PetscReal *)H, nc * nq * dim * dim * c);
274 }
275 PetscCheck(cell >= cStart && cell < cEnd, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %" PetscInt_FMT " not a cell [%" PetscInt_FMT ",%" PetscInt_FMT "), not implemented yet", cell, cStart, cEnd);
276 for (i = 0; i < nc * whol; i++) work[i] = dafield->cornerCoeffs[i];
277 for (j = 0; j < dim; j++) {
278 PetscReal e, d;
279 ijk[j] = (rem % cellsPer[j]);
280 rem /= cellsPer[j];
281
282 e = 2. * (ijk[j] + first[j] + 0.5) * stepPer[j] - 1.;
283 d = stepPer[j];
284 for (i = 0; i < half; i++) {
285 for (k = 0; k < nc; k++) {
286 cellCoeffs[i * nc + k] = work[2 * i * nc + k] * d;
287 cellCoeffs[(i + half) * nc + k] = work[2 * i * nc + k] * e + work[(2 * i + 1) * nc + k];
288 }
289 }
290 for (i = 0; i < whol * nc; i++) work[i] = cellCoeffs[i];
291 }
292 MultilinearEvaluate(dim, cellCoordRange, nc, cellCoeffs, dafield->work, nq, qs, datatype, cB, cD, cH);
293 }
294 if (!isStride) PetscCall(ISRestoreIndices(cellIS, &cells));
295 PetscCall(DMRestoreWorkArray(dm, (1 << dim) * nc, MPIU_SCALAR, &cellCoeffs));
296 #if defined(PETSC_USE_COMPLEX)
297 PetscCall(DMRestoreWorkArray(dm, nq * dim, MPIU_SCALAR, &qs));
298 #endif
299 PetscFunctionReturn(PETSC_SUCCESS);
300 }
301
DMFieldEvaluateFV_DA(DMField field,IS cellIS,PetscDataType datatype,void * B,void * D,void * H)302 static PetscErrorCode DMFieldEvaluateFV_DA(DMField field, IS cellIS, PetscDataType datatype, void *B, void *D, void *H)
303 {
304 PetscInt c, i, dim, cellsPer[3] = {0}, first[3] = {0};
305 PetscReal stepPer[3] = {0.};
306 DM dm;
307 DMDALocalInfo info;
308 PetscInt cStart, cEnd, numCells;
309 PetscInt nc;
310 PetscScalar *points;
311 DMField_DA *dafield;
312 PetscBool isStride;
313 const PetscInt *cells = NULL;
314 PetscInt sfirst = -1, stride = -1;
315
316 PetscFunctionBegin;
317 dafield = (DMField_DA *)field->data;
318 dm = field->dm;
319 nc = field->numComponents;
320 PetscCall(DMDAGetLocalInfo(dm, &info));
321 dim = info.dim;
322 stepPer[0] = 1. / info.mx;
323 stepPer[1] = 1. / info.my;
324 stepPer[2] = 1. / info.mz;
325 first[0] = info.gxs;
326 first[1] = info.gys;
327 first[2] = info.gzs;
328 cellsPer[0] = info.gxm;
329 cellsPer[1] = info.gym;
330 cellsPer[2] = info.gzm;
331 PetscCall(DMDAGetHeightStratum(dm, 0, &cStart, &cEnd));
332 PetscCall(ISGetLocalSize(cellIS, &numCells));
333 PetscCall(DMGetWorkArray(dm, dim * numCells, MPIU_SCALAR, &points));
334 PetscCall(PetscObjectTypeCompare((PetscObject)cellIS, ISSTRIDE, &isStride));
335 if (isStride) PetscCall(ISStrideGetInfo(cellIS, &sfirst, &stride));
336 else PetscCall(ISGetIndices(cellIS, &cells));
337 for (c = 0; c < numCells; c++) {
338 PetscInt cell = isStride ? (sfirst + c * stride) : cells[c];
339 PetscInt rem = cell;
340 PetscInt ijk[3] = {0};
341
342 PetscCheck(cell >= cStart && cell < cEnd, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %" PetscInt_FMT " not a cell [%" PetscInt_FMT ",%" PetscInt_FMT "), not implemented yet", cell, cStart, cEnd);
343 for (i = 0; i < dim; i++) {
344 ijk[i] = (rem % cellsPer[i]);
345 rem /= cellsPer[i];
346 points[dim * c + i] = (ijk[i] + first[i] + 0.5) * stepPer[i];
347 }
348 }
349 if (!isStride) PetscCall(ISRestoreIndices(cellIS, &cells));
350 MultilinearEvaluate(dim, dafield->coordRange, nc, dafield->cornerCoeffs, dafield->work, numCells, points, datatype, B, D, H);
351 PetscCall(DMRestoreWorkArray(dm, dim * numCells, MPIU_SCALAR, &points));
352 PetscFunctionReturn(PETSC_SUCCESS);
353 }
354
DMFieldGetDegree_DA(DMField field,IS pointIS,PetscInt * minDegree,PetscInt * maxDegree)355 static PetscErrorCode DMFieldGetDegree_DA(DMField field, IS pointIS, PetscInt *minDegree, PetscInt *maxDegree)
356 {
357 DM dm;
358 PetscInt dim, h, imin;
359
360 PetscFunctionBegin;
361 dm = field->dm;
362 PetscCall(ISGetMinMax(pointIS, &imin, NULL));
363 PetscCall(DMGetDimension(dm, &dim));
364 for (h = 0; h <= dim; h++) {
365 PetscInt hEnd;
366
367 PetscCall(DMDAGetHeightStratum(dm, h, NULL, &hEnd));
368 if (imin < hEnd) break;
369 }
370 dim -= h;
371 if (minDegree) *minDegree = 1;
372 if (maxDegree) *maxDegree = dim;
373 PetscFunctionReturn(PETSC_SUCCESS);
374 }
375
DMFieldCreateDefaultQuadrature_DA(DMField field,IS cellIS,PetscQuadrature * quad)376 static PetscErrorCode DMFieldCreateDefaultQuadrature_DA(DMField field, IS cellIS, PetscQuadrature *quad)
377 {
378 PetscInt h, dim, imax, imin;
379 DM dm;
380
381 PetscFunctionBegin;
382 dm = field->dm;
383 PetscCall(ISGetMinMax(cellIS, &imin, &imax));
384 PetscCall(DMGetDimension(dm, &dim));
385 *quad = NULL;
386 for (h = 0; h <= dim; h++) {
387 PetscInt hStart, hEnd;
388
389 PetscCall(DMDAGetHeightStratum(dm, h, &hStart, &hEnd));
390 if (imin >= hStart && imax < hEnd) break;
391 }
392 dim -= h;
393 if (dim > 0) PetscCall(PetscDTGaussTensorQuadrature(dim, 1, 1, -1.0, 1.0, quad));
394 PetscFunctionReturn(PETSC_SUCCESS);
395 }
396
DMFieldInitialize_DA(DMField field)397 static PetscErrorCode DMFieldInitialize_DA(DMField field)
398 {
399 DM dm;
400 Vec coords = NULL;
401 PetscInt dim, i, j, k;
402 DMField_DA *dafield = (DMField_DA *)field->data;
403
404 PetscFunctionBegin;
405 field->ops->destroy = DMFieldDestroy_DA;
406 field->ops->evaluate = DMFieldEvaluate_DA;
407 field->ops->evaluateFE = DMFieldEvaluateFE_DA;
408 field->ops->evaluateFV = DMFieldEvaluateFV_DA;
409 field->ops->getDegree = DMFieldGetDegree_DA;
410 field->ops->createDefaultQuadrature = DMFieldCreateDefaultQuadrature_DA;
411 field->ops->view = DMFieldView_DA;
412 dm = field->dm;
413 PetscCall(DMGetDimension(dm, &dim));
414 PetscCall(DMGetCoordinates(dm, &coords));
415 if (coords) {
416 PetscInt n;
417 const PetscScalar *array;
418 PetscReal mins[3][2] = {
419 {PETSC_MAX_REAL, PETSC_MAX_REAL},
420 {PETSC_MAX_REAL, PETSC_MAX_REAL},
421 {PETSC_MAX_REAL, PETSC_MAX_REAL}
422 };
423
424 PetscCall(VecGetLocalSize(coords, &n));
425 n /= dim;
426 PetscCall(VecGetArrayRead(coords, &array));
427 for (i = 0, k = 0; i < n; i++) {
428 for (j = 0; j < dim; j++, k++) {
429 PetscReal val = PetscRealPart(array[k]);
430
431 mins[j][0] = PetscMin(mins[j][0], val);
432 mins[j][1] = PetscMin(mins[j][1], -val);
433 }
434 }
435 PetscCall(VecRestoreArrayRead(coords, &array));
436 PetscCallMPI(MPIU_Allreduce(mins, &dafield->coordRange[0][0], 2 * dim, MPIU_REAL, MPI_MIN, PetscObjectComm((PetscObject)dm)));
437 for (j = 0; j < dim; j++) dafield->coordRange[j][1] = -dafield->coordRange[j][1];
438 } else {
439 for (j = 0; j < dim; j++) {
440 dafield->coordRange[j][0] = 0.;
441 dafield->coordRange[j][1] = 1.;
442 }
443 }
444 for (j = 0; j < dim; j++) {
445 PetscReal avg = 0.5 * (dafield->coordRange[j][1] + dafield->coordRange[j][0]);
446 PetscReal dif = 0.5 * (dafield->coordRange[j][1] - dafield->coordRange[j][0]);
447
448 dafield->coordRange[j][0] = avg;
449 dafield->coordRange[j][1] = dif;
450 }
451 PetscFunctionReturn(PETSC_SUCCESS);
452 }
453
DMFieldCreate_DA(DMField field)454 PETSC_INTERN PetscErrorCode DMFieldCreate_DA(DMField field)
455 {
456 DMField_DA *dafield;
457
458 PetscFunctionBegin;
459 PetscCall(PetscNew(&dafield));
460 field->data = dafield;
461 PetscCall(DMFieldInitialize_DA(field));
462 PetscFunctionReturn(PETSC_SUCCESS);
463 }
464
DMFieldCreateDA(DM dm,PetscInt nc,const PetscScalar * cornerValues,DMField * field)465 PetscErrorCode DMFieldCreateDA(DM dm, PetscInt nc, const PetscScalar *cornerValues, DMField *field)
466 {
467 DMField b;
468 DMField_DA *dafield;
469 PetscInt dim, nv, i, j, k;
470 PetscInt half;
471 PetscScalar *cv, *cf, *work;
472
473 PetscFunctionBegin;
474 PetscCall(DMFieldCreate(dm, nc, DMFIELD_VERTEX, &b));
475 PetscCall(DMFieldSetType(b, DMFIELDDA));
476 dafield = (DMField_DA *)b->data;
477 PetscCall(DMGetDimension(dm, &dim));
478 nv = (1 << dim) * nc;
479 PetscCall(PetscMalloc3(nv, &cv, nv, &cf, nv, &work));
480 for (i = 0; i < nv; i++) cv[i] = cornerValues[i];
481 for (i = 0; i < nv; i++) cf[i] = cv[i];
482 dafield->cornerVals = cv;
483 dafield->cornerCoeffs = cf;
484 dafield->work = work;
485 half = (1 << (dim - 1));
486 for (i = 0; i < dim; i++) {
487 PetscScalar *w;
488
489 w = work;
490 for (j = 0; j < half; j++) {
491 for (k = 0; k < nc; k++) w[j * nc + k] = 0.5 * (cf[(2 * j + 1) * nc + k] - cf[(2 * j) * nc + k]);
492 }
493 w = &work[j * nc];
494 for (j = 0; j < half; j++) {
495 for (k = 0; k < nc; k++) w[j * nc + k] = 0.5 * (cf[(2 * j) * nc + k] + cf[(2 * j + 1) * nc + k]);
496 }
497 for (j = 0; j < nv; j++) cf[j] = work[j];
498 }
499 *field = b;
500 PetscFunctionReturn(PETSC_SUCCESS);
501 }
502