Actual source code: cartesian.c
petsc-3.4.5 2014-06-29
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 = 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;
98: PetscBool iascii, isbinary, isdraw;
102: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
103: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERBINARY, &isbinary);
104: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERDRAW, &isdraw);
106: DMCartesianGetMesh(dm, m);
107: if (iascii) {
108: DMView_Cartesian_Ascii(m, viewer);
109: }
110: return(0);
111: }
115: PetscErrorCode DMCreateInterpolation_Cartesian(DM fineMesh, DM coarseMesh, Mat *interpolation, Vec *scaling)
116: {
117: ALE::Obj<ALE::CartesianMesh> coarse;
118: ALE::Obj<ALE::CartesianMesh> fine;
119: Mat P;
120: PetscErrorCode ierr;
123: DMCartesianGetMesh(fineMesh, fine);
124: DMCartesianGetMesh(coarseMesh, coarse);
125: #if 0
126: const ALE::Obj<ALE::Mesh::real_section_type>& coarseCoordinates = coarse->getRealSection("coordinates");
127: const ALE::Obj<ALE::Mesh::real_section_type>& fineCoordinates = fine->getRealSection("coordinates");
128: const ALE::Obj<ALE::Mesh::label_sequence>& vertices = fine->depthStratum(0);
129: const ALE::Obj<ALE::Mesh::real_section_type>& sCoarse = coarse->getRealSection("default");
130: const ALE::Obj<ALE::Mesh::real_section_type>& sFine = fine->getRealSection("default");
131: const ALE::Obj<ALE::Mesh::order_type>& coarseOrder = coarse->getFactory()->getGlobalOrder(coarse, "default", sCoarse);
132: const ALE::Obj<ALE::Mesh::order_type>& fineOrder = fine->getFactory()->getGlobalOrder(fine, "default", sFine);
133: const int dim = coarse->getDimension();
134: double *v0, *J, *invJ, detJ, *refCoords, *values;
135: #endif
137: MatCreate(fine->comm(), &P);
138: #if 0
139: MatSetSizes(P, sFine->size(), sCoarse->size(), PETSC_DETERMINE, PETSC_DETERMINE);
140: MatSetFromOptions(P);
141: PetscMalloc5(dim,double,&v0,dim*dim,double,&J,dim*dim,double,&invJ,dim,double,&refCoords,dim+1,double,&values);
142: for (ALE::Mesh::label_sequence::iterator v_iter = vertices->begin(); v_iter != vertices->end(); ++v_iter) {
143: const ALE::Mesh::real_section_type::value_type *coords = fineCoordinates->restrictPoint(*v_iter);
144: const ALE::Mesh::point_type coarseCell = coarse->locatePoint(coords);
146: coarse->computeElementGeometry(coarseCoordinates, coarseCell, v0, J, invJ, detJ);
147: for (int d = 0; d < dim; ++d) {
148: refCoords[d] = 0.0;
149: for (int e = 0; e < dim; ++e) {
150: refCoords[d] += invJ[d*dim+e]*(coords[e] - v0[e]);
151: }
152: refCoords[d] -= 1.0;
153: }
154: values[0] = 1.0/3.0 - (refCoords[0] + refCoords[1])/3.0;
155: values[1] = 0.5*(refCoords[0] + 1.0);
156: values[2] = 0.5*(refCoords[1] + 1.0);
157: updateOperatorGeneral(P, fine, sFine, fineOrder, *v_iter, sCoarse, coarseOrder, coarseCell, values, INSERT_VALUES);
158: }
159: PetscFree5(v0,J,invJ,refCoords,values);
160: #endif
161: MatAssemblyBegin(P, MAT_FINAL_ASSEMBLY);
162: MatAssemblyEnd(P, MAT_FINAL_ASSEMBLY);
163: *interpolation = P;
164: return(0);
165: }
169: PetscErrorCode DMRefine_Cartesian(DM mesh, MPI_Comm comm, DM *refinedMesh)
170: {
171: ALE::Obj<ALE::CartesianMesh> oldMesh;
172: PetscErrorCode ierr;
175: DMCartesianGetMesh(mesh, oldMesh);
176: DMCartesianCreate(comm, refinedMesh);
177: #if 0
178: ALE::Obj<ALE::Mesh> newMesh = ALE::Generator::refineMesh(oldMesh, refinementLimit, false);
179: MeshCartesianSetMesh(*refinedMesh, newMesh);
180: const ALE::Obj<ALE::CartesianMesh::real_section_type>& s = newMesh->getRealSection("default");
182: newMesh->setDiscretization(oldMesh->getDiscretization());
183: newMesh->setBoundaryCondition(oldMesh->getBoundaryCondition());
184: newMesh->setupField(s);
185: #else
186: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP, "Not yet implemented");
187: #endif
188: return(0);
189: }
193: PetscErrorCode DMCoarsen_Cartesian(DM mesh, MPI_Comm comm, DM *coarseMesh)
194: {
195: ALE::Obj<ALE::CartesianMesh> oldMesh;
196: PetscErrorCode ierr;
199: if (comm == MPI_COMM_NULL) {
200: PetscObjectGetComm((PetscObject)mesh,&comm);
201: }
202: DMCartesianGetMesh(mesh, oldMesh);
203: DMCartesianCreate(comm, coarseMesh);
204: #if 0
205: ALE::Obj<ALE::Mesh> newMesh = ALE::Generator::refineMesh(oldMesh, refinementLimit, false);
206: MeshCartesianSetMesh(*coarseMesh, newMesh);
207: const ALE::Obj<ALE::CartesianMesh::real_section_type>& s = newMesh->getRealSection("default");
209: newMesh->setDiscretization(oldMesh->getDiscretization());
210: newMesh->setBoundaryCondition(oldMesh->getBoundaryCondition());
211: newMesh->setupField(s);
212: #else
213: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP, "Not yet implemented");
214: #endif
215: return(0);
216: }
220: PetscErrorCode DMCartesianGetSectionReal(DM dm, const char name[], SectionReal *section)
221: {
222: ALE::Obj<ALE::CartesianMesh> m;
223: PetscErrorCode ierr;
226: DMCartesianGetMesh(dm, m);
227: SectionRealCreate(m->comm(), section);
228: PetscObjectSetName((PetscObject) *section, name);
229: #if 0
230: SectionRealSetSection(*section, m->getRealSection(std::string(name)));
231: SectionRealSetBundle(*section, m);
232: #endif
233: return(0);
234: }
238: PetscErrorCode DMSetFromOptions_Cartesian(DM dm)
239: {
244: PetscOptionsHead("DMCartesian Options");
245: /* Handle DMCartesian refinement */
246: /* Handle associated vectors */
247: PetscOptionsTail();
248: return(0);
249: }
253: PETSC_EXTERN 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>(NULL);
265: DMSetVecType(dm,VECSTANDARD);
267: dm->ops->globaltolocalbegin = 0;
268: dm->ops->globaltolocalend = 0;
269: dm->ops->localtoglobalbegin = 0;
270: dm->ops->localtoglobalend = 0;
271: dm->ops->createglobalvector = 0; /* DMCreateGlobalVector_Cartesian; */
272: dm->ops->createlocalvector = 0; /* DMCreateLocalVector_Cartesian; */
273: dm->ops->createinterpolation = DMCreateInterpolation_Cartesian;
274: dm->ops->getcoloring = 0;
275: dm->ops->creatematrix = 0; /* DMCreateMatrix_Cartesian; */
276: dm->ops->refine = DMRefine_Cartesian;
277: dm->ops->coarsen = DMCoarsen_Cartesian;
278: dm->ops->refinehierarchy = 0;
279: dm->ops->coarsenhierarchy = 0;
280: dm->ops->getinjection = 0;
281: dm->ops->getaggregates = 0;
282: dm->ops->destroy = DMDestroy_Cartesian;
283: dm->ops->view = DMView_Cartesian;
284: dm->ops->setfromoptions = DMSetFromOptions_Cartesian;
285: dm->ops->setup = 0;
286: return(0);
287: }
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: }