1 #include <petscdm.h> 2 #include <petscdmda.h> 3 #include <petscdmswarm.h> 4 #include <petsc/private/dmswarmimpl.h> 5 #include "../src/dm/impls/swarm/data_bucket.h" 6 7 PetscErrorCode private_PetscViewerCreate_XDMF(MPI_Comm comm, const char filename[], PetscViewer *v) { 8 long int *bytes; 9 PetscContainer container; 10 PetscViewer viewer; 11 12 PetscFunctionBegin; 13 PetscCall(PetscViewerCreate(comm, &viewer)); 14 PetscCall(PetscViewerSetType(viewer, PETSCVIEWERASCII)); 15 PetscCall(PetscViewerFileSetMode(viewer, FILE_MODE_WRITE)); 16 PetscCall(PetscViewerFileSetName(viewer, filename)); 17 18 PetscCall(PetscMalloc1(1, &bytes)); 19 bytes[0] = 0; 20 PetscCall(PetscContainerCreate(comm, &container)); 21 PetscCall(PetscContainerSetPointer(container, (void *)bytes)); 22 PetscCall(PetscObjectCompose((PetscObject)viewer, "XDMFViewerContext", (PetscObject)container)); 23 24 /* write xdmf header */ 25 PetscCall(PetscViewerASCIIPrintf(viewer, "<?xml version=\"1.0\" encoding=\"utf-8\"?>\n")); 26 PetscCall(PetscViewerASCIIPrintf(viewer, "<Xdmf xmlns:xi=\"http://www.w3.org/2001/XInclude/\" Version=\"2.99\">\n")); 27 /* write xdmf domain */ 28 PetscCall(PetscViewerASCIIPushTab(viewer)); 29 PetscCall(PetscViewerASCIIPrintf(viewer, "<Domain>\n")); 30 *v = viewer; 31 PetscFunctionReturn(0); 32 } 33 34 PetscErrorCode private_PetscViewerDestroy_XDMF(PetscViewer *v) { 35 PetscViewer viewer; 36 DM dm = NULL; 37 long int *bytes; 38 PetscContainer container = NULL; 39 40 PetscFunctionBegin; 41 if (!v) PetscFunctionReturn(0); 42 viewer = *v; 43 44 PetscCall(PetscObjectQuery((PetscObject)viewer, "DMSwarm", (PetscObject *)&dm)); 45 if (dm) { 46 PetscCall(PetscViewerASCIIPrintf(viewer, "</Grid>\n")); 47 PetscCall(PetscViewerASCIIPopTab(viewer)); 48 } 49 50 /* close xdmf header */ 51 PetscCall(PetscViewerASCIIPrintf(viewer, "</Domain>\n")); 52 PetscCall(PetscViewerASCIIPopTab(viewer)); 53 PetscCall(PetscViewerASCIIPrintf(viewer, "</Xdmf>\n")); 54 55 PetscCall(PetscObjectQuery((PetscObject)viewer, "XDMFViewerContext", (PetscObject *)&container)); 56 if (container) { 57 PetscCall(PetscContainerGetPointer(container, (void **)&bytes)); 58 PetscCall(PetscFree(bytes)); 59 PetscCall(PetscContainerDestroy(&container)); 60 } 61 PetscCall(PetscViewerDestroy(&viewer)); 62 *v = NULL; 63 PetscFunctionReturn(0); 64 } 65 66 PetscErrorCode private_CreateDataFileNameXDMF(const char filename[], char dfilename[]) { 67 char *ext; 68 PetscBool flg; 69 70 PetscFunctionBegin; 71 PetscCall(PetscStrrchr(filename, '.', &ext)); 72 PetscCall(PetscStrcmp("xmf", ext, &flg)); 73 if (flg) { 74 size_t len; 75 char viewername_minus_ext[PETSC_MAX_PATH_LEN]; 76 77 PetscCall(PetscStrlen(filename, &len)); 78 PetscCall(PetscStrncpy(viewername_minus_ext, filename, len - 2)); 79 PetscCall(PetscSNPrintf(dfilename, PETSC_MAX_PATH_LEN - 1, "%s_swarm_fields.pbin", viewername_minus_ext)); 80 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "File extension must by .xmf"); 81 PetscFunctionReturn(0); 82 } 83 84 PetscErrorCode private_DMSwarmView_XDMF(DM dm, PetscViewer viewer) { 85 PetscBool isswarm = PETSC_FALSE; 86 const char *viewername; 87 char datafile[PETSC_MAX_PATH_LEN]; 88 char *datafilename; 89 PetscViewer fviewer; 90 PetscInt k, ng, dim; 91 Vec dvec; 92 long int *bytes = NULL; 93 PetscContainer container = NULL; 94 const char *dmname; 95 96 PetscFunctionBegin; 97 PetscCall(PetscObjectQuery((PetscObject)viewer, "XDMFViewerContext", (PetscObject *)&container)); 98 if (container) { 99 PetscCall(PetscContainerGetPointer(container, (void **)&bytes)); 100 } else SETERRQ(PetscObjectComm((PetscObject)viewer), PETSC_ERR_SUP, "Valid to find attached data XDMFViewerContext"); 101 102 PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMSWARM, &isswarm)); 103 PetscCheck(isswarm, PetscObjectComm((PetscObject)viewer), PETSC_ERR_SUP, "Only valid for DMSwarm"); 104 105 PetscCall(PetscObjectCompose((PetscObject)viewer, "DMSwarm", (PetscObject)dm)); 106 107 PetscCall(PetscViewerASCIIPushTab(viewer)); 108 PetscCall(PetscObjectGetName((PetscObject)dm, &dmname)); 109 if (!dmname) { PetscCall(DMGetOptionsPrefix(dm, &dmname)); } 110 if (!dmname) { 111 PetscCall(PetscViewerASCIIPrintf(viewer, "<Grid Name=\"DMSwarm\" GridType=\"Uniform\">\n")); 112 } else { 113 PetscCall(PetscViewerASCIIPrintf(viewer, "<Grid Name=\"DMSwarm[%s]\" GridType=\"Uniform\">\n", dmname)); 114 } 115 116 /* create a sub-viewer for topology, geometry and all data fields */ 117 /* name is viewer.name + "_swarm_fields.pbin" */ 118 PetscCall(PetscViewerCreate(PetscObjectComm((PetscObject)viewer), &fviewer)); 119 PetscCall(PetscViewerSetType(fviewer, PETSCVIEWERBINARY)); 120 PetscCall(PetscViewerBinarySetSkipHeader(fviewer, PETSC_TRUE)); 121 PetscCall(PetscViewerBinarySetSkipInfo(fviewer, PETSC_TRUE)); 122 PetscCall(PetscViewerFileSetMode(fviewer, FILE_MODE_WRITE)); 123 124 PetscCall(PetscViewerFileGetName(viewer, &viewername)); 125 PetscCall(private_CreateDataFileNameXDMF(viewername, datafile)); 126 PetscCall(PetscViewerFileSetName(fviewer, datafile)); 127 PetscCall(PetscStrrchr(datafile, '/', &datafilename)); 128 129 PetscCall(DMSwarmGetSize(dm, &ng)); 130 131 /* write topology header */ 132 PetscCall(PetscViewerASCIIPushTab(viewer)); 133 PetscCall(PetscViewerASCIIPrintf(viewer, "<Topology Dimensions=\"%" PetscInt_FMT "\" TopologyType=\"Mixed\">\n", ng)); 134 PetscCall(PetscViewerASCIIPushTab(viewer)); 135 PetscCall(PetscViewerASCIIPrintf(viewer, "<DataItem Format=\"Binary\" Endian=\"Big\" DataType=\"Int\" Dimensions=\"%" PetscInt_FMT "\" Seek=\"%ld\">\n", ng * 3, bytes[0])); 136 PetscCall(PetscViewerASCIIPushTab(viewer)); 137 PetscCall(PetscViewerASCIIPrintf(viewer, "%s\n", datafilename)); 138 PetscCall(PetscViewerASCIIPopTab(viewer)); 139 PetscCall(PetscViewerASCIIPrintf(viewer, "</DataItem>\n")); 140 PetscCall(PetscViewerASCIIPopTab(viewer)); 141 PetscCall(PetscViewerASCIIPrintf(viewer, "</Topology>\n")); 142 PetscCall(PetscViewerASCIIPopTab(viewer)); 143 144 /* write topology data */ 145 for (k = 0; k < ng; k++) { 146 PetscInt pvertex[3]; 147 148 pvertex[0] = 1; 149 pvertex[1] = 1; 150 pvertex[2] = k; 151 PetscCall(PetscViewerBinaryWrite(fviewer, pvertex, 3, PETSC_INT)); 152 } 153 bytes[0] += sizeof(PetscInt) * ng * 3; 154 155 /* write geometry header */ 156 PetscCall(PetscViewerASCIIPushTab(viewer)); 157 PetscCall(DMGetDimension(dm, &dim)); 158 switch (dim) { 159 case 1: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "No support for 1D"); 160 case 2: PetscCall(PetscViewerASCIIPrintf(viewer, "<Geometry Type=\"XY\">\n")); break; 161 case 3: PetscCall(PetscViewerASCIIPrintf(viewer, "<Geometry Type=\"XYZ\">\n")); break; 162 } 163 PetscCall(PetscViewerASCIIPushTab(viewer)); 164 PetscCall(PetscViewerASCIIPrintf(viewer, "<DataItem Format=\"Binary\" Endian=\"Big\" DataType=\"Float\" Precision=\"8\" Dimensions=\"%" PetscInt_FMT " %" PetscInt_FMT "\" Seek=\"%ld\">\n", ng, dim, bytes[0])); 165 PetscCall(PetscViewerASCIIPushTab(viewer)); 166 PetscCall(PetscViewerASCIIPrintf(viewer, "%s\n", datafilename)); 167 PetscCall(PetscViewerASCIIPopTab(viewer)); 168 PetscCall(PetscViewerASCIIPrintf(viewer, "</DataItem>\n")); 169 PetscCall(PetscViewerASCIIPopTab(viewer)); 170 PetscCall(PetscViewerASCIIPrintf(viewer, "</Geometry>\n")); 171 PetscCall(PetscViewerASCIIPopTab(viewer)); 172 173 /* write geometry data */ 174 PetscCall(DMSwarmCreateGlobalVectorFromField(dm, DMSwarmPICField_coor, &dvec)); 175 PetscCall(VecView(dvec, fviewer)); 176 PetscCall(DMSwarmDestroyGlobalVectorFromField(dm, DMSwarmPICField_coor, &dvec)); 177 bytes[0] += sizeof(PetscReal) * ng * dim; 178 179 PetscCall(PetscViewerDestroy(&fviewer)); 180 PetscFunctionReturn(0); 181 } 182 183 PetscErrorCode private_VecView_Swarm_XDMF(Vec x, PetscViewer viewer) { 184 long int *bytes = NULL; 185 PetscContainer container = NULL; 186 const char *viewername; 187 char datafile[PETSC_MAX_PATH_LEN]; 188 char *datafilename; 189 PetscViewer fviewer; 190 PetscInt N, bs; 191 const char *vecname; 192 char fieldname[PETSC_MAX_PATH_LEN]; 193 194 PetscFunctionBegin; 195 PetscCall(PetscObjectQuery((PetscObject)viewer, "XDMFViewerContext", (PetscObject *)&container)); 196 PetscCheck(container, PetscObjectComm((PetscObject)viewer), PETSC_ERR_SUP, "Unable to find attached data XDMFViewerContext"); 197 PetscCall(PetscContainerGetPointer(container, (void **)&bytes)); 198 PetscCall(PetscViewerFileGetName(viewer, &viewername)); 199 PetscCall(private_CreateDataFileNameXDMF(viewername, datafile)); 200 201 /* re-open a sub-viewer for all data fields */ 202 /* name is viewer.name + "_swarm_fields.pbin" */ 203 PetscCall(PetscViewerCreate(PetscObjectComm((PetscObject)viewer), &fviewer)); 204 PetscCall(PetscViewerSetType(fviewer, PETSCVIEWERBINARY)); 205 PetscCall(PetscViewerBinarySetSkipHeader(fviewer, PETSC_TRUE)); 206 PetscCall(PetscViewerBinarySetSkipInfo(fviewer, PETSC_TRUE)); 207 PetscCall(PetscViewerFileSetMode(fviewer, FILE_MODE_APPEND)); 208 PetscCall(PetscViewerFileSetName(fviewer, datafile)); 209 PetscCall(PetscStrrchr(datafile, '/', &datafilename)); 210 211 PetscCall(VecGetSize(x, &N)); 212 PetscCall(VecGetBlockSize(x, &bs)); 213 N = N / bs; 214 PetscCall(PetscObjectGetName((PetscObject)x, &vecname)); 215 if (!vecname) { 216 PetscCall(PetscSNPrintf(fieldname, PETSC_MAX_PATH_LEN - 1, "swarmfield_%d", ((PetscObject)x)->tag)); 217 } else { 218 PetscCall(PetscSNPrintf(fieldname, PETSC_MAX_PATH_LEN - 1, "%s", vecname)); 219 } 220 221 /* write data header */ 222 PetscCall(PetscViewerASCIIPushTab(viewer)); 223 PetscCall(PetscViewerASCIIPrintf(viewer, "<Attribute Center=\"Node\" Name=\"%s\" Type=\"None\">\n", fieldname)); 224 PetscCall(PetscViewerASCIIPushTab(viewer)); 225 if (bs == 1) { 226 PetscCall(PetscViewerASCIIPrintf(viewer, "<DataItem Format=\"Binary\" Endian=\"Big\" DataType=\"Float\" Precision=\"8\" Dimensions=\"%" PetscInt_FMT "\" Seek=\"%ld\">\n", N, bytes[0])); 227 } else { 228 PetscCall(PetscViewerASCIIPrintf(viewer, "<DataItem Format=\"Binary\" Endian=\"Big\" DataType=\"Float\" Precision=\"8\" Dimensions=\"%" PetscInt_FMT " %" PetscInt_FMT "\" Seek=\"%ld\">\n", N, bs, bytes[0])); 229 } 230 PetscCall(PetscViewerASCIIPushTab(viewer)); 231 PetscCall(PetscViewerASCIIPrintf(viewer, "%s\n", datafilename)); 232 PetscCall(PetscViewerASCIIPopTab(viewer)); 233 PetscCall(PetscViewerASCIIPrintf(viewer, "</DataItem>\n")); 234 PetscCall(PetscViewerASCIIPopTab(viewer)); 235 PetscCall(PetscViewerASCIIPrintf(viewer, "</Attribute>\n")); 236 PetscCall(PetscViewerASCIIPopTab(viewer)); 237 238 /* write data */ 239 PetscCall(VecView(x, fviewer)); 240 bytes[0] += sizeof(PetscReal) * N * bs; 241 242 PetscCall(PetscViewerDestroy(&fviewer)); 243 PetscFunctionReturn(0); 244 } 245 246 PetscErrorCode private_ISView_Swarm_XDMF(IS is, PetscViewer viewer) { 247 long int *bytes = NULL; 248 PetscContainer container = NULL; 249 const char *viewername; 250 char datafile[PETSC_MAX_PATH_LEN]; 251 char *datafilename; 252 PetscViewer fviewer; 253 PetscInt N, bs; 254 const char *vecname; 255 char fieldname[PETSC_MAX_PATH_LEN]; 256 257 PetscFunctionBegin; 258 PetscCall(PetscObjectQuery((PetscObject)viewer, "XDMFViewerContext", (PetscObject *)&container)); 259 PetscCheck(container, PetscObjectComm((PetscObject)viewer), PETSC_ERR_SUP, "Unable to find attached data XDMFViewerContext"); 260 PetscCall(PetscContainerGetPointer(container, (void **)&bytes)); 261 PetscCall(PetscViewerFileGetName(viewer, &viewername)); 262 PetscCall(private_CreateDataFileNameXDMF(viewername, datafile)); 263 264 /* re-open a sub-viewer for all data fields */ 265 /* name is viewer.name + "_swarm_fields.pbin" */ 266 PetscCall(PetscViewerCreate(PetscObjectComm((PetscObject)viewer), &fviewer)); 267 PetscCall(PetscViewerSetType(fviewer, PETSCVIEWERBINARY)); 268 PetscCall(PetscViewerBinarySetSkipHeader(fviewer, PETSC_TRUE)); 269 PetscCall(PetscViewerBinarySetSkipInfo(fviewer, PETSC_TRUE)); 270 PetscCall(PetscViewerFileSetMode(fviewer, FILE_MODE_APPEND)); 271 PetscCall(PetscViewerFileSetName(fviewer, datafile)); 272 PetscCall(PetscStrrchr(datafile, '/', &datafilename)); 273 274 PetscCall(ISGetSize(is, &N)); 275 PetscCall(ISGetBlockSize(is, &bs)); 276 N = N / bs; 277 PetscCall(PetscObjectGetName((PetscObject)is, &vecname)); 278 if (!vecname) { 279 PetscCall(PetscSNPrintf(fieldname, PETSC_MAX_PATH_LEN - 1, "swarmfield_%d", ((PetscObject)is)->tag)); 280 } else { 281 PetscCall(PetscSNPrintf(fieldname, PETSC_MAX_PATH_LEN - 1, "%s", vecname)); 282 } 283 284 /* write data header */ 285 PetscCall(PetscViewerASCIIPushTab(viewer)); 286 PetscCall(PetscViewerASCIIPrintf(viewer, "<Attribute Center=\"Node\" Name=\"%s\" Type=\"None\">\n", fieldname)); 287 PetscCall(PetscViewerASCIIPushTab(viewer)); 288 if (bs == 1) { 289 PetscCall(PetscViewerASCIIPrintf(viewer, "<DataItem Format=\"Binary\" Endian=\"Big\" DataType=\"Int\" Precision=\"4\" Dimensions=\"%" PetscInt_FMT "\" Seek=\"%ld\">\n", N, bytes[0])); 290 } else { 291 PetscCall(PetscViewerASCIIPrintf(viewer, "<DataItem Format=\"Binary\" Endian=\"Big\" DataType=\"Int\" Precision=\"4\" Dimensions=\"%" PetscInt_FMT " %" PetscInt_FMT "\" Seek=\"%ld\">\n", N, bs, bytes[0])); 292 } 293 PetscCall(PetscViewerASCIIPushTab(viewer)); 294 PetscCall(PetscViewerASCIIPrintf(viewer, "%s\n", datafilename)); 295 PetscCall(PetscViewerASCIIPopTab(viewer)); 296 PetscCall(PetscViewerASCIIPrintf(viewer, "</DataItem>\n")); 297 PetscCall(PetscViewerASCIIPopTab(viewer)); 298 PetscCall(PetscViewerASCIIPrintf(viewer, "</Attribute>\n")); 299 PetscCall(PetscViewerASCIIPopTab(viewer)); 300 301 /* write data */ 302 PetscCall(ISView(is, fviewer)); 303 bytes[0] += sizeof(PetscInt) * N * bs; 304 305 PetscCall(PetscViewerDestroy(&fviewer)); 306 PetscFunctionReturn(0); 307 } 308 309 /*@C 310 DMSwarmViewFieldsXDMF - Write a selection of DMSwarm fields to an XDMF3 file 311 312 Collective on dm 313 314 Input parameters: 315 + dm - the DMSwarm 316 . filename - the file name of the XDMF file (must have the extension .xmf) 317 . nfields - the number of fields to write into the XDMF file 318 - field_name_list - array of length nfields containing the textual name of fields to write 319 320 Level: beginner 321 322 Notes: 323 Only fields registered with data type PETSC_DOUBLE or PETSC_INT can be written into the file 324 325 .seealso: `DMSwarmViewXDMF()` 326 @*/ 327 PETSC_EXTERN PetscErrorCode DMSwarmViewFieldsXDMF(DM dm, const char filename[], PetscInt nfields, const char *field_name_list[]) { 328 Vec dvec; 329 PetscInt f, N; 330 PetscViewer viewer; 331 332 PetscFunctionBegin; 333 PetscCall(private_PetscViewerCreate_XDMF(PetscObjectComm((PetscObject)dm), filename, &viewer)); 334 PetscCall(private_DMSwarmView_XDMF(dm, viewer)); 335 PetscCall(DMSwarmGetLocalSize(dm, &N)); 336 for (f = 0; f < nfields; f++) { 337 void *data; 338 PetscDataType type; 339 340 PetscCall(DMSwarmGetField(dm, field_name_list[f], NULL, &type, &data)); 341 PetscCall(DMSwarmRestoreField(dm, field_name_list[f], NULL, &type, &data)); 342 if (type == PETSC_DOUBLE) { 343 PetscCall(DMSwarmCreateGlobalVectorFromField(dm, field_name_list[f], &dvec)); 344 PetscCall(PetscObjectSetName((PetscObject)dvec, field_name_list[f])); 345 PetscCall(private_VecView_Swarm_XDMF(dvec, viewer)); 346 PetscCall(DMSwarmDestroyGlobalVectorFromField(dm, field_name_list[f], &dvec)); 347 } else if (type == PETSC_INT) { 348 IS is; 349 const PetscInt *idx; 350 351 PetscCall(DMSwarmGetField(dm, field_name_list[f], NULL, &type, &data)); 352 idx = (const PetscInt *)data; 353 354 PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)dm), N, idx, PETSC_USE_POINTER, &is)); 355 PetscCall(PetscObjectSetName((PetscObject)is, field_name_list[f])); 356 PetscCall(private_ISView_Swarm_XDMF(is, viewer)); 357 PetscCall(ISDestroy(&is)); 358 PetscCall(DMSwarmRestoreField(dm, field_name_list[f], NULL, &type, &data)); 359 } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Can only write PETSC_INT and PETSC_DOUBLE"); 360 } 361 PetscCall(private_PetscViewerDestroy_XDMF(&viewer)); 362 PetscFunctionReturn(0); 363 } 364 365 /*@C 366 DMSwarmViewXDMF - Write DMSwarm fields to an XDMF3 file 367 368 Collective on dm 369 370 Input parameters: 371 + dm - the DMSwarm 372 - filename - the file name of the XDMF file (must have the extension .xmf) 373 374 Level: beginner 375 376 Notes: 377 Only fields user registered with data type PETSC_DOUBLE or PETSC_INT will be written into the file 378 379 Developer Notes: 380 This should be removed and replaced with the standard use of PetscViewer 381 382 .seealso: `DMSwarmViewFieldsXDMF()` 383 @*/ 384 PETSC_EXTERN PetscErrorCode DMSwarmViewXDMF(DM dm, const char filename[]) { 385 DM_Swarm *swarm = (DM_Swarm *)dm->data; 386 Vec dvec; 387 PetscInt f; 388 PetscViewer viewer; 389 390 PetscFunctionBegin; 391 PetscCall(private_PetscViewerCreate_XDMF(PetscObjectComm((PetscObject)dm), filename, &viewer)); 392 PetscCall(private_DMSwarmView_XDMF(dm, viewer)); 393 for (f = 4; f < swarm->db->nfields; f++) { /* only examine user defined fields - the first 4 are internally created by DMSwarmPIC */ 394 DMSwarmDataField field; 395 396 /* query field type - accept all those of type PETSC_DOUBLE */ 397 field = swarm->db->field[f]; 398 if (field->petsc_type == PETSC_DOUBLE) { 399 PetscCall(DMSwarmCreateGlobalVectorFromField(dm, field->name, &dvec)); 400 PetscCall(PetscObjectSetName((PetscObject)dvec, field->name)); 401 PetscCall(private_VecView_Swarm_XDMF(dvec, viewer)); 402 PetscCall(DMSwarmDestroyGlobalVectorFromField(dm, field->name, &dvec)); 403 } else if (field->petsc_type == PETSC_INT) { 404 IS is; 405 PetscInt N; 406 const PetscInt *idx; 407 void *data; 408 409 PetscCall(DMSwarmGetLocalSize(dm, &N)); 410 PetscCall(DMSwarmGetField(dm, field->name, NULL, NULL, &data)); 411 idx = (const PetscInt *)data; 412 413 PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)dm), N, idx, PETSC_USE_POINTER, &is)); 414 PetscCall(PetscObjectSetName((PetscObject)is, field->name)); 415 PetscCall(private_ISView_Swarm_XDMF(is, viewer)); 416 PetscCall(ISDestroy(&is)); 417 PetscCall(DMSwarmRestoreField(dm, field->name, NULL, NULL, &data)); 418 } 419 } 420 PetscCall(private_PetscViewerDestroy_XDMF(&viewer)); 421 PetscFunctionReturn(0); 422 } 423