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 ? 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
