## Re: [Libmesh-devel] Vector-valued FE

 Re: [Libmesh-devel] Vector-valued FE From: Roy Stogner - 2012-05-24 15:51:30 ```On Thu, 10 May 2012, Derek Gaston wrote: > As long as we're on this subject let me throw a few things out there.... Between illness and work and the renumbering discussion I let this slip... > What I've run into though is that there is no good way to get gradient > evaluation vectorized because we are working with vectors of Gradient > objects... meaning that the values themselves are not contiguous in memory. >  So when we're doing something like this: > >   for (unsigned int i=0; i < num_dofs; i++) >   { >     soln_local = current_solution(dof_indices[i]); > >     for (unsigned int qp=0; qp < nqp; qp++) >       grad_u_qp->add_scaled(grad_phi[i][qp], soln_local); >   } > > that inner for loop can't be vectorized because those grad_phi[i][qp] > objects aren't contiguous memory. I've seen this problem with non-gradients: we want qp to be incremented in the outer loop, because there's solution variables to be assembled and/or nonlinear constitutive terms to be evaluated at each quadrature point, but then the inner loop is not vectorizable because we're iterating over phi[i][qp] rather than phi[qp][i]. I've got a couple ideas that might fix that (and give you your flat vectors, now that you mention it) in the long run. But... since when are grad_phi_i[qp] objects not continuous memory? If the underlying data there isn't an n_qp*LIBMESH_DIM array of Real, then something's going wrong under the hood. > We are now supporting multiply dimensioned pieces of mesh in the same mesh > (ie we have 1D elements and 3D elements in the same mesh... but not sharing > nodes (although they do transfer information back and forth)).  In order to > do that seamlessly we have to copy the shape function values in and out of > internal data structures.  It would be awesome if we could _pass in_ the > vectors to be filled by libMesh instead of having libMesh own the data and > pass us back references.  Then we could avoid a lot of copying by passing > the same vector to be filled into multiple FE objects (that are of different > dimension) and asking them to fill those vectors when we're on the > corresponding elements... Hmmm... I'm still not seeing how this is useful (why make a copy rather than just pass around the original?) but it wouldn't be too hard to support. Probably falls underneath "vectorize the shape function calls in reinit()" on my eternally growing TODO list, though. --- Roy```