Actual source code: cartesian.c
petsc-3.3-p7 2013-05-11
1: #include <petsc-private/meshimpl.h> /*I "petscdmmesh.h" I*/
2: #include <CartesianSieve.hh>
6: /*@C
7: DMCartesianGetMesh - Gets the internal mesh object
9: Not collective
11: Input Parameter:
12: . dm - the mesh object
14: Output Parameter:
15: . m - the internal mesh object
17: Level: advanced
19: .seealso MeshCreate(), MeshCartesianSetMesh()
20: @*/
21: PetscErrorCode DMCartesianGetMesh(DM dm, ALE::Obj<ALE::CartesianMesh>& m)
22: {
23: DM_Cartesian *c = (DM_Cartesian *) dm->data;
27: m = c->m;
28: return(0);
29: }
33: /*@C
34: DMCartesianSetMesh - Sets the internal mesh object
36: Not collective
38: Input Parameters:
39: + mesh - the mesh object
40: - m - the internal mesh object
42: Level: advanced
44: .seealso MeshCreate(), MeshCartesianGetMesh()
45: @*/
46: PetscErrorCode DMCartesianSetMesh(DM dm, const ALE::Obj<ALE::CartesianMesh>& m)
47: {
48: DM_Cartesian *c = (DM_Cartesian *) dm->data;
52: c->m = m;
53: return(0);
54: }
58: PetscErrorCode DMDestroy_Cartesian(DM dm)
59: {
60: DM_Cartesian *c = (DM_Cartesian *) dm->data;
63: c->m = PETSC_NULL;
64: return(0);
65: }
69: PetscErrorCode DMView_Cartesian_Ascii(const ALE::Obj<ALE::CartesianMesh>& mesh, PetscViewer viewer)
70: {
71: PetscViewerFormat format;
72: PetscErrorCode ierr;
75: PetscViewerGetFormat(viewer, &format);
76: if (format == PETSC_VIEWER_ASCII_VTK) {
77: #if 0
78: VTKViewer::writeHeader(viewer);
79: VTKViewer::writeVertices(mesh, viewer);
80: VTKViewer::writeElements(mesh, viewer);
81: #endif
82: } else {
83: int dim = mesh->getDimension();
85: PetscViewerASCIIPrintf(viewer, "Mesh in %d dimensions:\n", dim);
86: /* FIX: Need to globalize */
87: PetscViewerASCIIPrintf(viewer, " %d vertices\n", mesh->getSieve()->getNumVertices());
88: PetscViewerASCIIPrintf(viewer, " %d cells\n", mesh->getSieve()->getNumCells());
89: }
90: PetscViewerFlush(viewer);
91: return(0);
92: }
94: PetscErrorCode DMView_Cartesian(DM dm, PetscViewer viewer)
95: {
96: ALE::Obj<ALE::CartesianMesh> m;
97: PetscBool iascii, isbinary, isdraw;
101: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
102: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERBINARY, &isbinary);
103: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERDRAW, &isdraw);
105: DMCartesianGetMesh(dm, m);
106: if (iascii){
107: DMView_Cartesian_Ascii(m, viewer);
108: } else if (isbinary) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP, "Binary viewer not implemented for Cartesian Mesh");
109: else if (isdraw) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP, "Draw viewer not implemented for Cartesian Mesh");
110: else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Viewer type %s not supported by this mesh object", ((PetscObject)viewer)->type_name);
111: return(0);
112: }
116: PetscErrorCode DMCreateInterpolation_Cartesian(DM fineMesh, DM coarseMesh, Mat *interpolation, Vec *scaling)
117: {
118: ALE::Obj<ALE::CartesianMesh> coarse;
119: ALE::Obj<ALE::CartesianMesh> fine;
120: Mat P;
121: PetscErrorCode ierr;
124: DMCartesianGetMesh(fineMesh, fine);
125: DMCartesianGetMesh(coarseMesh, coarse);
126: #if 0
127: const ALE::Obj<ALE::Mesh::real_section_type>& coarseCoordinates = coarse->getRealSection("coordinates");
128: const ALE::Obj<ALE::Mesh::real_section_type>& fineCoordinates = fine->getRealSection("coordinates");
129: const ALE::Obj<ALE::Mesh::label_sequence>& vertices = fine->depthStratum(0);
130: const ALE::Obj<ALE::Mesh::real_section_type>& sCoarse = coarse->getRealSection("default");
131: const ALE::Obj<ALE::Mesh::real_section_type>& sFine = fine->getRealSection("default");
132: const ALE::Obj<ALE::Mesh::order_type>& coarseOrder = coarse->getFactory()->getGlobalOrder(coarse, "default", sCoarse);
133: const ALE::Obj<ALE::Mesh::order_type>& fineOrder = fine->getFactory()->getGlobalOrder(fine, "default", sFine);
134: const int dim = coarse->getDimension();
135: double *v0, *J, *invJ, detJ, *refCoords, *values;
136: #endif
138: MatCreate(fine->comm(), &P);
139: #if 0
140: MatSetSizes(P, sFine->size(), sCoarse->size(), PETSC_DETERMINE, PETSC_DETERMINE);
141: MatSetFromOptions(P);
142: PetscMalloc5(dim,double,&v0,dim*dim,double,&J,dim*dim,double,&invJ,dim,double,&refCoords,dim+1,double,&values);
143: for(ALE::Mesh::label_sequence::iterator v_iter = vertices->begin(); v_iter != vertices->end(); ++v_iter) {
144: const ALE::Mesh::real_section_type::value_type *coords = fineCoordinates->restrictPoint(*v_iter);
145: const ALE::Mesh::point_type coarseCell = coarse->locatePoint(coords);
147: coarse->computeElementGeometry(coarseCoordinates, coarseCell, v0, J, invJ, detJ);
148: for(int d = 0; d < dim; ++d) {
149: refCoords[d] = 0.0;
150: for(int e = 0; e < dim; ++e) {
151: refCoords[d] += invJ[d*dim+e]*(coords[e] - v0[e]);
152: }
153: refCoords[d] -= 1.0;
154: }
155: values[0] = 1.0/3.0 - (refCoords[0] + refCoords[1])/3.0;
156: values[1] = 0.5*(refCoords[0] + 1.0);
157: values[2] = 0.5*(refCoords[1] + 1.0);
158: updateOperatorGeneral(P, fine, sFine, fineOrder, *v_iter, sCoarse, coarseOrder, coarseCell, values, INSERT_VALUES);
159: }
160: PetscFree5(v0,J,invJ,refCoords,values);
161: #endif
162: MatAssemblyBegin(P, MAT_FINAL_ASSEMBLY);
163: MatAssemblyEnd(P, MAT_FINAL_ASSEMBLY);
164: *interpolation = P;
165: return(0);
166: }
170: PetscErrorCode DMRefine_Cartesian(DM mesh, MPI_Comm comm, DM *refinedMesh)
171: {
172: ALE::Obj<ALE::CartesianMesh> oldMesh;
173: PetscErrorCode ierr;
176: DMCartesianGetMesh(mesh, oldMesh);
177: DMCartesianCreate(comm, refinedMesh);
178: #if 0
179: ALE::Obj<ALE::Mesh> newMesh = ALE::Generator::refineMesh(oldMesh, refinementLimit, false);
180: MeshCartesianSetMesh(*refinedMesh, newMesh);
181: const ALE::Obj<ALE::CartesianMesh::real_section_type>& s = newMesh->getRealSection("default");
183: newMesh->setDiscretization(oldMesh->getDiscretization());
184: newMesh->setBoundaryCondition(oldMesh->getBoundaryCondition());
185: newMesh->setupField(s);
186: #else
187: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP, "Not yet implemented");
188: #endif
189: return(0);
190: }
194: PetscErrorCode DMCoarsen_Cartesian(DM mesh, MPI_Comm comm, DM *coarseMesh)
195: {
196: ALE::Obj<ALE::CartesianMesh> oldMesh;
197: PetscErrorCode ierr;
200: if (comm == MPI_COMM_NULL) comm = ((PetscObject)mesh)->comm;
201: DMCartesianGetMesh(mesh, oldMesh);
202: DMCartesianCreate(comm, coarseMesh);
203: #if 0
204: ALE::Obj<ALE::Mesh> newMesh = ALE::Generator::refineMesh(oldMesh, refinementLimit, false);
205: MeshCartesianSetMesh(*coarseMesh, newMesh);
206: const ALE::Obj<ALE::CartesianMesh::real_section_type>& s = newMesh->getRealSection("default");
208: newMesh->setDiscretization(oldMesh->getDiscretization());
209: newMesh->setBoundaryCondition(oldMesh->getBoundaryCondition());
210: newMesh->setupField(s);
211: #else
212: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP, "Not yet implemented");
213: #endif
214: return(0);
215: }
219: PetscErrorCode DMCartesianGetSectionReal(DM dm, const char name[], SectionReal *section)
220: {
221: ALE::Obj<ALE::CartesianMesh> m;
222: PetscErrorCode ierr;
225: DMCartesianGetMesh(dm, m);
226: SectionRealCreate(m->comm(), section);
227: PetscObjectSetName((PetscObject) *section, name);
228: #if 0
229: SectionRealSetSection(*section, m->getRealSection(std::string(name)));
230: SectionRealSetBundle(*section, m);
231: #endif
232: return(0);
233: }
237: PetscErrorCode DMSetFromOptions_Cartesian(DM dm)
238: {
243: PetscOptionsHead("DMCartesian Options");
244: /* Handle DMCartesian refinement */
245: /* Handle associated vectors */
246: PetscOptionsTail();
247: return(0);
248: }
250: EXTERN_C_BEGIN
253: PetscErrorCode DMCreate_Cartesian(DM dm)
254: {
255: DM_Cartesian *mesh;
260: PetscNewLog(dm, DM_Cartesian, &mesh);
261: dm->data = mesh;
263: new(&mesh->m) ALE::Obj<ALE::CartesianMesh>(PETSC_NULL);
265: PetscStrallocpy(VECSTANDARD, &dm->vectype);
266: dm->ops->globaltolocalbegin = 0;
267: dm->ops->globaltolocalend = 0;
268: dm->ops->localtoglobalbegin = 0;
269: dm->ops->localtoglobalend = 0;
270: dm->ops->createglobalvector = 0; /* DMCreateGlobalVector_Cartesian; */
271: dm->ops->createlocalvector = 0; /* DMCreateLocalVector_Cartesian; */
272: dm->ops->createinterpolation = DMCreateInterpolation_Cartesian;
273: dm->ops->getcoloring = 0;
274: dm->ops->creatematrix = 0; /* DMCreateMatrix_Cartesian; */
275: dm->ops->refine = DMRefine_Cartesian;
276: dm->ops->coarsen = DMCoarsen_Cartesian;
277: dm->ops->refinehierarchy = 0;
278: dm->ops->coarsenhierarchy = 0;
279: dm->ops->getinjection = 0;
280: dm->ops->getaggregates = 0;
281: dm->ops->destroy = DMDestroy_Cartesian;
282: dm->ops->view = DMView_Cartesian;
283: dm->ops->setfromoptions = DMSetFromOptions_Cartesian;
284: dm->ops->setup = 0;
285: return(0);
286: }
287: EXTERN_C_END
291: /*@
292: DMCartesianCreate - Creates a DMCartesian object.
294: Collective on MPI_Comm
296: Input Parameter:
297: . comm - The communicator for the DMCartesian object
299: Output Parameter:
300: . mesh - The DMCartesian object
302: Level: beginner
304: .keywords: DMCartesian, create
305: @*/
306: PetscErrorCode DMCartesianCreate(MPI_Comm comm, DM *mesh)
307: {
312: DMCreate(comm, mesh);
313: DMSetType(*mesh, DMCARTESIAN);
314: return(0);
315: }