On Fri, Jul 10, 2015 at 2:56 PM, Marco Zocca <[email protected]> wrote:
> Thank you for the reply; > > let me improve my question: in e.g. a FEM code, we need: > > > *) a mesh table :: [element -> [face -> [edge -> [node] ] ] ] > > *) one or more reference elements with associated basis over elements > and/or faces and Jacobians, or a quadrature rule to represent > integration over the real-space elements, for each element > > *) one or more scalar/tensor coefficients per element > > > Assembly (in the simplest case) of the global matrix is then a single > loop over elements. > > > I would like to use the highest-level primitives for this task (I am > tackling a stochastic PDE-constrained optimization problem, so > complexity control and modularity is paramount ), could you provide > some hints ? > I try to organize things this way in SNES ex12, but there is a lot of associated infrastructure, so it might not be all that easy to read. However, I orthogonalize mesh topology from data layout from element specification from form integration from assembly from solvers, so that you can change them independently. Matt > Thank you again and kind regards, > Marco > > > > > > > >> I'm looking for some non-optimized PETSc code; namely, I struggle a > >> bit with generalizing the provided examples. On one hand inlined > >> elementary operations help keep track of the FLOPs but make for > >> hard-to read code. > >> > > > > > You can use MatSetValuesStencil() as is done here: > > > > > http://www.mcs.anl.gov/petsc/petsc-current/src/ksp/ksp/examples/tutorials/ex50.c.html > > > > on line 162. > > > >> > > >> , removed comments for space) we are interested in forming a matrix A > >> from the values contained in a vector x; i.e. some 'map' operation > >> that sweeps x and prepares rows of A at a time according to some fixed > >> stencil, with separate treatment of boundary values . > >> > >> Isn't there an idiom to abstract out this functionality? It's one of > >> the most common operations in numerical PDE codes: > > > > > > Actually, this is only common is toy example. Almost no real PDEs have a > > fixed stencil > > because you have material coefficients, or nonlinearities, or > constraints, > > etc. > > > > Matt > -- 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
