On Fri, Aug 1, 2014 at 6:01 PM, Justin Chang <[email protected]> wrote:
> Never mind, I sort of figured it out. What I was kind of looking for was > within the DMPlexMatSetClosure(...) function (plex.c). Specifically, these > lines and the private functions within: > > 5383: DMGetWorkArray > <http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/DM/DMGetWorkArray.html#DMGetWorkArray>(dm, > numIndices, PETSC_INT, &indices);5384: if (numFields) {5385: for (p = > 0; p < numPoints*2; p += 2) {5386: PetscInt > <http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Sys/PetscInt.html#PetscInt> > o = points[p+1];5387: PetscSectionGetOffset > <http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/IS/PetscSectionGetOffset.html#PetscSectionGetOffset>(globalSection, > points[p], &globalOff);5388: indicesPointFields_private(section, > points[p], globalOff < 0 ? -(globalOff+1) : globalOff, offsets, PETSC_FALSE > <http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Sys/PETSC_FALSE.html#PETSC_FALSE>, > o, indices);5389: }5390: } else {5391: for (p = 0, off = 0; p < > numPoints*2; p += 2) {5392: PetscInt > <http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Sys/PetscInt.html#PetscInt> > o = points[p+1];5393: PetscSectionGetOffset > <http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/IS/PetscSectionGetOffset.html#PetscSectionGetOffset>(globalSection, > points[p], &globalOff);5394: indicesPoint_private(section, points[p], > globalOff < 0 ? -(globalOff+1) : globalOff, &off, PETSC_FALSE > <http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Sys/PETSC_FALSE.html#PETSC_FALSE>, > o, indices);5395: }5396: } > > If I understand this part of the code correctly, it returns the global > indices based on the global offsets; if a certain global dof happens to be > constrained the indices will be negative. All I want is the ability to > extract global column indices because each row of my equality constraint > matrix will be element specific. Is there an existing function that does > something like this, or do I have to write it myself? > > I think you have to write that one yourself. The closest function is probably PetscErrorCode DMPlexMatGetClosureIndicesRefined(DM dmf, PetscSection fsection, PetscSection globalFSection, DM dmc, PetscSection csection, PetscSection globalCSection, PetscInt point, PetscInt cindices[], PetscInt findices[]) Matt > On Fri, Aug 1, 2014 at 1:40 PM, Matthew Knepley <[email protected]> wrote: > >> On Fri, Aug 1, 2014 at 2:38 PM, Justin Chang <[email protected]> wrote: >> >>> As far as I know this equality constraint should exist outside of the >>> large square matrix. >>> >>> What I was thinking of doing is manually creating an *Nele* by *Ndof* >>> matrix, each row will have the same number of non-zeros. Since the size of >>> *Ndof* is simply the global number of free dofs formulated by the >>> PetscSection, I want to know what's the best way to map the free dofs into >>> the global column indices [0 *Ndof*) >>> >>> For example, in src/dm/impls/plex/examples/tutorials/ex1.c where a box >>> mesh with three fields is created, there are a total of 41 dofs to which 8 >>> of them are constrained. If no constraints were set, then the size of u >>> (i.e. *Ndof*) would be 41, but with constraints set it is 33. When you >>> view the constrained global vector u, it lists all the dofs in the order as >>> specified by the PetscSection but skips all the constrained ones. Would I >>> have to use something like DMSet/GetDefaultGlobalSection to do the mapping? >>> >> >> I do not understand what you are asking. However, you can just use the >> PetscSection for the Ndof field to give you >> the local sizes for your matrix. >> >> Matt >> >> >>> Thanks, >>> Justin >>> >>> >>> On Fri, Aug 1, 2014 at 7:26 AM, Matthew Knepley <[email protected]> >>> wrote: >>> >>>> On Wed, Jul 30, 2014 at 5:01 PM, Justin Chang <[email protected]> >>>> wrote: >>>> >>>>> That's good to know thanks. I have another question somewhat related. >>>>> >>>>> In that quadratic programming problem I have shown earlier, my K >>>>> matrix is of global size *Ndof* by *Ndof* where *Ndof* is number of >>>>> free dofs (for instance, if I have a 2x2 quadrilateral mesh and want to >>>>> evaluate concentration at all the verticies, I will have nine total dofs >>>>> but if I have Dirichlet constraints all along the boundary, only my center >>>>> vertex is free so *Ndof* = 1. Thus when I simply call >>>>> DMCreateMatrix(...) it would give me a 1x1 matrix.) >>>>> >>>>> However, I want my M constraint matrix to be of size *Nele* by *Ndof* >>>>> where *Nele* is the total number of elements. How would I create this >>>>> matrix? Because ultimately, I am subjecting all of my free u's to >>>>> *Nele* number of equality constraints, and if possible I want this M >>>>> matrix to be compatible with the layout from Vec u. >>>>> >>>> >>>> This is not hard in principle. I used to have code that did this, but I >>>> took it out during the rewrite because no one >>>> was using it and it was complicated. It just generalizes the code in >>>> DMPlexPreallocateOperator() to take two >>>> PetscSections instead of one. >>>> >>>> Are you sure that this does not belong inside a large square matrix for >>>> the entire problem? >>>> >>>> Thanks, >>>> >>>> Matt >>>> >>>> >>>>> >>>>> On Wed, Jul 30, 2014 at 2:03 PM, Matthew Knepley <[email protected]> >>>>> wrote: >>>>> >>>>>> On Wed, Jul 30, 2014 at 2:53 PM, Justin Chang <[email protected]> >>>>>> wrote: >>>>>> >>>>>>> Hi, >>>>>>> >>>>>>> This might be a silly question, but I am developing a finite element >>>>>>> code to solve the non-negative diffusion equation. I am using the DMPlex >>>>>>> structure and plan on subjecting the problem to a bounded constraint >>>>>>> optimization to ensure non-negative concentrations. Say for instance, I >>>>>>> want to solve: >>>>>>> >>>>>>> min || K*u - F ||^2 >>>>>>> subject to u >= 0 >>>>>>> M*u = 0 >>>>>>> >>>>>>> Where K is the Jacobian, F the residual, u the concentration, and M >>>>>>> some equality constraint matrix. I know how to do this via Tao, but my >>>>>>> question is can Tao handle Mats and Vecs created via DMPlex? Because >>>>>>> from >>>>>>> the SNES and TS examples I have seen in ex12, ex62, and ex11, it seems >>>>>>> there are solvers tailored to handle DMPlex created data structures. >>>>>>> >>>>>> >>>>>> Yes, TAO can handle objects created by DMPlex since they are just Vec >>>>>> and Mat structures. The additions >>>>>> to ex12, ex62, and ex11 concern the callbacks from the solver to the >>>>>> constructions routines for residual and >>>>>> Jacobian. Right now, none of the callbacks are automatic for TAO so >>>>>> you would have to code them yourself, >>>>>> probably by just calling the routines from SNES or TS so it should >>>>>> not be that hard. We are working to get >>>>>> everything integrated as it is in the other solvers. >>>>>> >>>>>> Thanks, >>>>>> >>>>>> Matt >>>>>> >>>>>> >>>>>>> Thanks, >>>>>>> Justin >>>>>>> >>>>>> >>>>>> >>>>>> >>>>>> -- >>>>>> What most experimenters take for granted before they begin their >>>>>> experiments is infinitely more interesting than any results to which >>>>>> their >>>>>> experiments lead. >>>>>> -- Norbert Wiener >>>>>> >>>>> >>>>> >>>> >>>> >>>> -- >>>> What most experimenters take for granted before they begin their >>>> experiments is infinitely more interesting than any results to which their >>>> experiments lead. >>>> -- Norbert Wiener >>>> >>> >>> >> >> >> -- >> What most experimenters take for granted before they begin their >> experiments is infinitely more interesting than any results to which their >> experiments lead. >> -- Norbert Wiener >> > > -- What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead. -- Norbert Wiener
