Actual source code: dmmbmg.cxx

  1: #include <petsc/private/dmmbimpl.h>
  2: #include <petscdmmoab.h>
  3: #include <MBTagConventions.hpp>
  4: #include <moab/NestedRefine.hpp>

  6: /*@C
  7:   DMMoabGenerateHierarchy - Generate a multi-level uniform refinement hierarchy
  8:   by succesively refining a coarse mesh, already defined in the DM object
  9:   provided by the user.

 11:   Collective

 13:   Input Parameter:
 14: . dmb  - The DMMoab object

 16:   Output Parameter:
 17: + nlevels   - The number of levels of refinement needed to generate the hierarchy
 18: - ldegrees  - The degree of refinement at each level in the hierarchy

 20:   Level: beginner

 22: @*/
 23: PetscErrorCode DMMoabGenerateHierarchy(DM dm, PetscInt nlevels, PetscInt *ldegrees)
 24: {
 25:   DM_Moab        *dmmoab;
 26:   PetscErrorCode  ierr;
 27:   moab::ErrorCode merr;
 28:   PetscInt *pdegrees, ilevel;
 29:   std::vector<moab::EntityHandle> hsets;

 33:   dmmoab = (DM_Moab*)(dm)->data;

 35:   if (!ldegrees) {
 36:     PetscMalloc1(nlevels, &pdegrees);
 37:     for (ilevel = 0; ilevel < nlevels; ilevel++) pdegrees[ilevel] = 2; /* default = Degree 2 refinement */
 38:   }
 39:   else pdegrees = ldegrees;

 41:   /* initialize set level refinement data for hierarchy */
 42:   dmmoab->nhlevels = nlevels;

 44:   /* Instantiate the nested refinement class */
 45: #ifdef MOAB_HAVE_MPI
 46:   dmmoab->hierarchy = new moab::NestedRefine(dynamic_cast<moab::Core*>(dmmoab->mbiface), dmmoab->pcomm, dmmoab->fileset);
 47: #else
 48:   dmmoab->hierarchy = new moab::NestedRefine(dynamic_cast<moab::Core*>(dmmoab->mbiface), NULL, dmmoab->fileset);
 49: #endif

 51:   PetscMalloc1(nlevels + 1, &dmmoab->hsets);

 53:   /* generate the mesh hierarchy */
 54:   merr = dmmoab->hierarchy->generate_mesh_hierarchy(nlevels, pdegrees, hsets, false);MBERRNM(merr);

 56: #ifdef MOAB_HAVE_MPI
 57:   if (dmmoab->pcomm->size() > 1) {
 58:     merr = dmmoab->hierarchy->exchange_ghosts(hsets, dmmoab->nghostrings);MBERRNM(merr);
 59:   }
 60: #endif

 62:   /* copy the mesh sets for nested refinement hierarchy */
 63:   dmmoab->hsets[0] = hsets[0];
 64:   for (ilevel = 1; ilevel <= nlevels; ilevel++)
 65:   {
 66:     dmmoab->hsets[ilevel] = hsets[ilevel];

 68: #ifdef MOAB_HAVE_MPI
 69:     merr = dmmoab->pcomm->assign_global_ids(hsets[ilevel], dmmoab->dim, 0, false, true, false);MBERRNM(merr);
 70: #endif

 72:     /* Update material and other geometric tags from parent to child sets */
 73:     merr = dmmoab->hierarchy->update_special_tags(ilevel, hsets[ilevel]);MBERRNM(merr);
 74:   }

 76:   hsets.clear();
 77:   if (!ldegrees) {
 78:     PetscFree(pdegrees);
 79:   }
 80:   return(0);
 81: }

 83: /*@C
 84:   DMRefineHierarchy_Moab - Generate a multi-level DM hierarchy
 85:   by succesively refining a coarse mesh.

 87:   Collective

 89:   Input Parameter:
 90: . dm  - The DMMoab object

 92:   Output Parameter:
 93: + nlevels   - The number of levels of refinement needed to generate the hierarchy
 94: - dmf  - The DM objects after successive refinement of the hierarchy

 96:   Level: beginner

 98: @*/
 99: PETSC_EXTERN PetscErrorCode  DMRefineHierarchy_Moab(DM dm, PetscInt nlevels, DM dmf[])
100: {
101:   PetscErrorCode  ierr;
102:   PetscInt        i;


106:   DMRefine(dm, PetscObjectComm((PetscObject)dm), &dmf[0]);
107:   for (i = 1; i < nlevels; i++) {
108:     DMRefine(dmf[i - 1], PetscObjectComm((PetscObject)dm), &dmf[i]);
109:   }
110:   return(0);
111: }

113: /*@C
114:   DMCoarsenHierarchy_Moab - Generate a multi-level DM hierarchy
115:   by succesively coarsening a refined mesh.

117:   Collective

119:   Input Parameter:
120: . dm  - The DMMoab object

122:   Output Parameter:
123: + nlevels   - The number of levels of refinement needed to generate the hierarchy
124: - dmc  - The DM objects after successive coarsening of the hierarchy

126:   Level: beginner

128: @*/
129: PETSC_EXTERN PetscErrorCode DMCoarsenHierarchy_Moab(DM dm, PetscInt nlevels, DM dmc[])
130: {
131:   PetscErrorCode  ierr;
132:   PetscInt        i;


136:   DMCoarsen(dm, PetscObjectComm((PetscObject)dm), &dmc[0]);
137:   for (i = 1; i < nlevels; i++) {
138:     DMCoarsen(dmc[i - 1], PetscObjectComm((PetscObject)dm), &dmc[i]);
139:   }
140:   return(0);
141: }

143: PETSC_EXTERN PetscErrorCode DMMoab_Compute_NNZ_From_Connectivity(DM, PetscInt*, PetscInt*, PetscInt*, PetscInt*, PetscBool);

145: /*@C
146:   DMCreateInterpolation_Moab - Generate the interpolation operators to transform
147:   operators (matrices, vectors) from parent level to child level as defined by
148:   the DM inputs provided by the user.

150:   Collective

152:   Input Parameter:
153: + dm1  - The DMMoab object
154: - dm2  - the second, finer DMMoab object

156:   Output Parameter:
157: + interpl  - The interpolation operator for transferring data between the levels
158: - vec      - The scaling vector (optional)

160:   Level: developer

162: @*/
163: PETSC_EXTERN PetscErrorCode DMCreateInterpolation_Moab(DM dmp, DM dmc, Mat* interpl, Vec* vec)
164: {
165:   DM_Moab         *dmbp, *dmbc;
166:   PetscErrorCode   ierr;
167:   moab::ErrorCode  merr;
168:   PetscInt         dim;
169:   PetscReal        factor;
170:   PetscInt         innz, *nnz, ionz, *onz;
171:   PetscInt         nlsizp, nlsizc, nlghsizp, ngsizp, ngsizc;
172:   const PetscBool  use_consistent_bases=PETSC_TRUE;

177:   dmbp = (DM_Moab*)(dmp)->data;
178:   dmbc = (DM_Moab*)(dmc)->data;
179:   nlsizp = dmbp->nloc;// *dmb1->numFields;
180:   nlsizc = dmbc->nloc;// *dmb2->numFields;
181:   ngsizp = dmbp->n; // *dmb1->numFields;
182:   ngsizc = dmbc->n; // *dmb2->numFields;
183:   nlghsizp = (dmbp->nloc + dmbp->nghost); // *dmb1->numFields;

185:   // Columns = Parent DoFs ;  Rows = Child DoFs
186:   // Interpolation matrix: \sum_{i=1}^P Owned(Child) * (Owned(Parent) + Ghosted(Parent))
187:   // Size: nlsizc * nlghsizp
188:   PetscInfo4(NULL, "Creating interpolation matrix %D X %D to apply transformation between levels %D -> %D.\n", ngsizc, nlghsizp, dmbp->hlevel, dmbc->hlevel);

190:   DMGetDimension(dmp, &dim);

192:   /* allocate the nnz, onz arrays based on block size and local nodes */
193:   PetscCalloc2(nlsizc, &nnz, nlsizc, &onz);

195:   /* Loop through the local elements and compute the relation between the current parent and the refined_level. */
196:   for (moab::Range::iterator iter = dmbc->vowned->begin(); iter != dmbc->vowned->end(); iter++) {

198:     const moab::EntityHandle vhandle = *iter;
199:     /* define local variables */
200:     moab::EntityHandle parent;
201:     std::vector<moab::EntityHandle> adjs;
202:     moab::Range     found;

204:     /* store the vertex DoF number */
205:     const int ldof = dmbc->lidmap[vhandle - dmbc->seqstart];

207:     /* Get adjacency information for current vertex - i.e., all elements of dimension (dim) that connects
208:        to the current vertex. We can then decipher if a vertex is ghosted or not and compute the
209:        non-zero pattern accordingly. */
210:     merr = dmbc->hierarchy->get_adjacencies(vhandle, dmbc->dim, adjs);MBERRNM(merr);

212:     /* loop over vertices and update the number of connectivity */
213:     for (unsigned jter = 0; jter < adjs.size(); jter++) {

215:       const moab::EntityHandle jhandle = adjs[jter];

217:       /* Get the relation between the current (coarse) parent and its corresponding (finer) children elements */
218:       merr = dmbc->hierarchy->child_to_parent(jhandle, dmbc->hlevel, dmbp->hlevel, &parent);MBERRNM(merr);

220:       /* Get connectivity information in canonical ordering for the local element */
221:       std::vector<moab::EntityHandle> connp;
222:       merr = dmbp->hierarchy->get_connectivity(parent, dmbp->hlevel, connp);MBERRNM(merr);

224:       for (unsigned ic=0; ic < connp.size(); ++ic) {

226:         /* loop over each element connected to the adjacent vertex and update as needed */
227:         /* find the truly user-expected layer of ghosted entities to decipher NNZ pattern */
228:         if (found.find(connp[ic]) != found.end()) continue; /* make sure we don't double count shared vertices */
229:         if (dmbp->vghost->find(connp[ic]) != dmbp->vghost->end()) onz[ldof]++; /* update out-of-proc onz */
230:         else nnz[ldof]++; /* else local vertex */
231:         found.insert(connp[ic]);
232:       }
233:     }
234:   }

236:   for (int i = 0; i < nlsizc; i++)
237:     nnz[i] += 1; /* self count the node */

239:   ionz = onz[0];
240:   innz = nnz[0];
241:   for (int tc = 0; tc < nlsizc; tc++) {
242:     // check for maximum allowed sparsity = fully dense
243:     nnz[tc] = std::min(nlsizp, nnz[tc]);
244:     onz[tc] = std::min(ngsizp - nlsizp, onz[tc]);

246:     PetscInfo3(NULL, "  %d: NNZ = %d, ONZ = %d\n", tc, nnz[tc], onz[tc]);

248:     innz = (innz < nnz[tc] ? nnz[tc] : innz);
249:     ionz = (ionz < onz[tc] ? onz[tc] : ionz);
250:   }

252:   /* create interpolation matrix */
253:   MatCreate(PetscObjectComm((PetscObject)dmc), interpl);
254:   MatSetSizes(*interpl, nlsizc, nlsizp, ngsizc, ngsizp);
255:   MatSetType(*interpl, MATAIJ);
256:   MatSetFromOptions(*interpl);

258:   MatSeqAIJSetPreallocation(*interpl, innz, nnz);
259:   MatMPIAIJSetPreallocation(*interpl, innz, nnz, ionz, onz);

261:   /* clean up temporary memory */
262:   PetscFree2(nnz, onz);

264:   /* set up internal matrix data-structures */
265:   MatSetUp(*interpl);

267:   /* Define variables for assembly */
268:   std::vector<moab::EntityHandle> children;
269:   std::vector<moab::EntityHandle> connp, connc;
270:   std::vector<PetscReal> pcoords, ccoords, values_phi;

272:   if (use_consistent_bases) {
273:     const moab::EntityHandle ehandle = dmbp->elocal->front();

275:     merr = dmbp->hierarchy->parent_to_child(ehandle, dmbp->hlevel, dmbc->hlevel, children);MBERRNM(merr);

277:     /* Get connectivity and coordinates of the parent vertices */
278:     merr = dmbp->hierarchy->get_connectivity(ehandle, dmbp->hlevel, connp);MBERRNM(merr);
279:     merr = dmbc->mbiface->get_connectivity(&children[0], children.size(), connc);MBERRNM(merr);

281:     std::vector<PetscReal> natparam(3*connc.size(), 0.0);
282:     pcoords.resize(connp.size() * 3);
283:     ccoords.resize(connc.size() * 3);
284:     values_phi.resize(connp.size()*connc.size());
285:     /* Get coordinates for connectivity entities in canonical order for both coarse and finer levels */
286:     merr = dmbp->hierarchy->get_coordinates(&connp[0], connp.size(), dmbp->hlevel, &pcoords[0]);MBERRNM(merr);
287:     merr = dmbc->hierarchy->get_coordinates(&connc[0], connc.size(), dmbc->hlevel, &ccoords[0]);MBERRNM(merr);

289:     /* Set values: For each DOF in coarse grid cell, set the contribution or PHI evaluated at each fine grid DOF point */
290:     for (unsigned tc = 0; tc < connc.size(); tc++) {
291:       const PetscInt offset = tc * 3;

293:       /* Scale ccoords relative to pcoords */
294:       DMMoabPToRMapping(dim, connp.size(), &pcoords[0], &ccoords[offset], &natparam[offset], &values_phi[connp.size()*tc]);
295:     }
296:   }
297:   else {
298:     factor = std::pow(2.0 /*degree_P_for_refinement*/, (dmbc->hlevel - dmbp->hlevel) * dmbp->dim * 1.0);
299:   }

301:   /* TODO: Decipher the correct non-zero pattern. There is still some issue with onz allocation */
302:   MatSetOption(*interpl, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);

304:   /* Loop through the remaining vertices. These vertices appear only on the current refined_level. */
305:   for (moab::Range::iterator iter = dmbp->elocal->begin(); iter != dmbp->elocal->end(); iter++) {

307:     const moab::EntityHandle ehandle = *iter;

309:     /* Get the relation between the current (coarse) parent and its corresponding (finer) children elements */
310:     children.clear();
311:     connc.clear();
312:     merr = dmbp->hierarchy->parent_to_child(ehandle, dmbp->hlevel, dmbc->hlevel, children);MBERRNM(merr);

314:     /* Get connectivity and coordinates of the parent vertices */
315:     merr = dmbp->hierarchy->get_connectivity(ehandle, dmbp->hlevel, connp);MBERRNM(merr);
316:     merr = dmbc->mbiface->get_connectivity(&children[0], children.size(), connc);MBERRNM(merr);

318:     pcoords.resize(connp.size() * 3);
319:     ccoords.resize(connc.size() * 3);
320:     /* Get coordinates for connectivity entities in canonical order for both coarse and finer levels */
321:     merr = dmbp->hierarchy->get_coordinates(&connp[0], connp.size(), dmbp->hlevel, &pcoords[0]);MBERRNM(merr);
322:     merr = dmbc->hierarchy->get_coordinates(&connc[0], connc.size(), dmbc->hlevel, &ccoords[0]);MBERRNM(merr);

324:     std::vector<int> dofsp(connp.size()), dofsc(connc.size());
325:     /* TODO: specific to scalar system - use GetDofs */
326:     DMMoabGetDofsBlocked(dmp, connp.size(), &connp[0], &dofsp[0]);
327:     DMMoabGetDofsBlocked(dmc, connc.size(), &connc[0], &dofsc[0]);

329:     /* Compute the actual interpolation weights when projecting solution/residual between levels */
330:     if (use_consistent_bases) {

332:       /* Use the cached values of natural parameteric coordinates and basis pre-evaluated.
333:          We are making an assumption here that UMR used in GMG to generate the hierarchy uses
334:          the same template for all elements; This will fail for mixed element meshes (TRI/QUAD).

336:          TODO: Fix the above assumption by caching data for families (especially for Tets and mixed meshes)
337:       */

339:       /* Set values: For each DOF in coarse grid cell, set the contribution or PHI evaluated at each fine grid DOF point */
340:       for (unsigned tc = 0; tc < connc.size(); tc++) {
341:         /* TODO: Check if we should be using INSERT_VALUES instead */
342:         MatSetValues(*interpl, 1, &dofsc[tc], connp.size(), &dofsp[0], &values_phi[connp.size()*tc], ADD_VALUES);
343:       }
344:     }
345:     else {
346:       /* Compute the interpolation weights by determining distance of 1-ring
347:          neighbor vertices from current vertex

349:          This should be used only when FEM basis is not used for the discretization.
350:          Else, the consistent interface to compute the basis function for interpolation
351:          between the levels should be evaluated correctly to preserve convergence of GMG.
352:          Shephard's basis will be terrible for any unsmooth problems.
353:       */
354:       values_phi.resize(connp.size());
355:       for (unsigned tc = 0; tc < connc.size(); tc++) {

357:         PetscReal normsum = 0.0;
358:         for (unsigned tp = 0; tp < connp.size(); tp++) {
359:           values_phi[tp] = 0.0;
360:           for (unsigned k = 0; k < 3; k++)
361:             values_phi[tp] += std::pow(pcoords[tp * 3 + k] - ccoords[k + tc * 3], dim);
362:           if (values_phi[tp] < 1e-12) {
363:             values_phi[tp] = 1e12;
364:           }
365:           else {
366:             //values_phi[tp] = std::pow(values_phi[tp], -1.0/dim);
367:             values_phi[tp] = std::pow(values_phi[tp], -1.0);
368:             normsum += values_phi[tp];
369:           }
370:         }
371:         for (unsigned tp = 0; tp < connp.size(); tp++) {
372:           if (values_phi[tp] > 1e11)
373:             values_phi[tp] = factor * 0.5 / connp.size();
374:           else
375:             values_phi[tp] = factor * values_phi[tp] * 0.5 / (connp.size() * normsum);
376:         }
377:         MatSetValues(*interpl, 1, &dofsc[tc], connp.size(), &dofsp[0], &values_phi[0], ADD_VALUES);
378:       }
379:     }
380:   }
381:   if (vec) *vec = NULL;
382:   MatAssemblyBegin(*interpl, MAT_FINAL_ASSEMBLY);
383:   MatAssemblyEnd(*interpl, MAT_FINAL_ASSEMBLY);
384:   return(0);
385: }

387: /*@C
388:   DMCreateInjection_Moab - Generate a multi-level uniform refinement hierarchy
389:   by succesively refining a coarse mesh, already defined in the DM object
390:   provided by the user.

392:   Collective

394:   Input Parameter:
395: . dmb  - The DMMoab object

397:   Output Parameter:
398: + nlevels   - The number of levels of refinement needed to generate the hierarchy
399: - ldegrees  - The degree of refinement at each level in the hierarchy

401:   Level: beginner

403: @*/
404: PETSC_EXTERN PetscErrorCode DMCreateInjection_Moab(DM dm1, DM dm2, VecScatter* ctx)
405: {
406:   //DM_Moab        *dmmoab;

411:   //dmmoab = (DM_Moab*)(dm1)->data;

413:   PetscPrintf(PETSC_COMM_WORLD, "[DMCreateInjection_Moab] :: Placeholder\n");
414:   return(0);
415: }

417: static PetscErrorCode DMMoab_UMR_Private(DM dm, MPI_Comm comm, PetscBool refine, DM *dmref)
418: {
419:   PetscErrorCode  ierr;
420:   PetscInt        i, dim;
421:   DM              dm2;
422:   moab::ErrorCode merr;
423:   DM_Moab        *dmb = (DM_Moab*)dm->data, *dd2;


429:   if ((dmb->hlevel == dmb->nhlevels && refine) || (dmb->hlevel == 0 && !refine)) {
430:     if (dmb->hlevel + 1 > dmb->nhlevels && refine) PetscInfo2(NULL, "Invalid multigrid refinement hierarchy level specified (%D). MOAB UMR max levels = %D. Creating a NULL object.\n", dmb->hlevel + 1, dmb->nhlevels);
431:     if (dmb->hlevel - 1 < 0 && !refine) PetscInfo1(NULL, "Invalid multigrid coarsen hierarchy level specified (%D). Creating a NULL object.\n", dmb->hlevel - 1);
432:     *dmref = PETSC_NULL;
433:     return(0);
434:   }

436:   DMMoabCreate(PetscObjectComm((PetscObject)dm), &dm2);
437:   dd2 = (DM_Moab*)dm2->data;

439:   dd2->mbiface = dmb->mbiface;
440: #ifdef MOAB_HAVE_MPI
441:   dd2->pcomm = dmb->pcomm;
442: #endif
443:   dd2->icreatedinstance = PETSC_FALSE;
444:   dd2->nghostrings = dmb->nghostrings;

446:   /* set the new level based on refinement/coarsening */
447:   if (refine) {
448:     dd2->hlevel = dmb->hlevel + 1;
449:   }
450:   else {
451:     dd2->hlevel = dmb->hlevel - 1;
452:   }

454:   /* Copy the multilevel hierarchy pointers in MOAB */
455:   dd2->hierarchy = dmb->hierarchy;
456:   dd2->nhlevels = dmb->nhlevels;
457:   PetscMalloc1(dd2->nhlevels + 1, &dd2->hsets);
458:   for (i = 0; i <= dd2->nhlevels; i++) {
459:     dd2->hsets[i] = dmb->hsets[i];
460:   }
461:   dd2->fileset = dd2->hsets[dd2->hlevel];

463:   /* do the remaining initializations for DMMoab */
464:   dd2->bs = dmb->bs;
465:   dd2->numFields = dmb->numFields;
466:   dd2->rw_dbglevel = dmb->rw_dbglevel;
467:   dd2->partition_by_rank = dmb->partition_by_rank;
468:   PetscStrcpy(dd2->extra_read_options, dmb->extra_read_options);
469:   PetscStrcpy(dd2->extra_write_options, dmb->extra_write_options);
470:   dd2->read_mode = dmb->read_mode;
471:   dd2->write_mode = dmb->write_mode;

473:   /* set global ID tag handle */
474:   DMMoabSetLocalToGlobalTag(dm2, dmb->ltog_tag);

476:   merr = dd2->mbiface->tag_get_handle(MATERIAL_SET_TAG_NAME, dd2->material_tag);MBERRNM(merr);

478:   DMSetOptionsPrefix(dm2, ((PetscObject)dm)->prefix);
479:   DMGetDimension(dm, &dim);
480:   DMSetDimension(dm2, dim);

482:   /* allow overloaded (user replaced) operations to be inherited by refinement clones */
483:   dm2->ops->creatematrix = dm->ops->creatematrix;

485:   /* copy fill information if given */
486:   DMMoabSetBlockFills(dm2, dmb->dfill, dmb->ofill);

488:   /* copy vector type information */
489:   DMSetMatType(dm2, dm->mattype);
490:   DMSetVecType(dm2, dm->vectype);
491:   dd2->numFields = dmb->numFields;
492:   if (dmb->numFields) {
493:     DMMoabSetFieldNames(dm2, dmb->numFields, dmb->fieldNames);
494:   }

496:   DMSetFromOptions(dm2);

498:   /* recreate Dof numbering for the refined DM and make sure the distribution is correctly populated */
499:   DMSetUp(dm2);

501:   *dmref = dm2;
502:   return(0);
503: }


506: /*@C
507:   DMRefine_Moab - Generate a multi-level uniform refinement hierarchy
508:   by succesively refining a coarse mesh, already defined in the DM object
509:   provided by the user.

511:   Collective on dm

513:   Input Parameter:
514: + dm  - The DMMoab object
515: - comm - the communicator to contain the new DM object (or MPI_COMM_NULL)

517:   Output Parameter:
518: . dmf - the refined DM, or NULL

520:   Note: If no refinement was done, the return value is NULL

522:   Level: developer

524: @*/
525: PETSC_EXTERN PetscErrorCode DMRefine_Moab(DM dm, MPI_Comm comm, DM* dmf)
526: {
527:   PetscErrorCode  ierr;


532:   DMMoab_UMR_Private(dm, comm, PETSC_TRUE, dmf);
533:   return(0);
534: }

536: /*@C
537:   DMCoarsen_Moab - Generate a multi-level uniform refinement hierarchy
538:   by succesively refining a coarse mesh, already defined in the DM object
539:   provided by the user.

541:   Collective on dm

543:   Input Parameter:
544: + dm  - The DMMoab object
545: - comm - the communicator to contain the new DM object (or MPI_COMM_NULL)

547:   Output Parameter:
548: . dmf - the coarsened DM, or NULL

550:   Note: If no coarsening was done, the return value is NULL

552:   Level: developer

554: @*/
555: PETSC_EXTERN PetscErrorCode DMCoarsen_Moab(DM dm, MPI_Comm comm, DM* dmc)
556: {
557:   PetscErrorCode  ierr;


562:   DMMoab_UMR_Private(dm, comm, PETSC_FALSE, dmc);
563:   return(0);
564: }