Dear James, Andriy is right of course, the gradient are computed in Getfem computing the analytical gradient of the shape functions on each element. The function getfem::compute_gradient interpolate the gradient on a given (Lagrange) finite element method which have to be a non-continuous one since the gradient are not continuous from an element to another.
Extracting the linear operator corresponding to the function getfem::compute_gradient as a matrix is not impossible of course. This means a matrix whose columns correspond to the derivative of each shape function. For the moment, such a function is not available in Getfem but it would correspond to something very similar to getfem::compute_gradient. Yves. Le 13/05/2012 14:42, Andriy Andreykiv a écrit : > Dear James, > > I'm not 100% sure in this specific case, but I would expect that > Getfem computes gradients not by finite differences, but by using the > analytical gradients of the > shape functions. Doing per element something like this: Grad U = Grad > Shape Func X U. This is a very common practice in FE if I'm not mistaken. > > Best regards, > Andriy > > On 13 May 2012 14:26, James Zhou <[email protected] > <mailto:[email protected]>> wrote: > > Thank you Andriy. The linear operator I meant is global. The > code you provided is very clear, but I am still wondering how the > gradient of displacement is computed? In your code, it > corresponds to "getfem::compute_gradient(mf_u, mf_strain, U, GRAD);". > > I have looked the method compute_gradient in > src/getfem/getfem_derivatives.h, but didn't quite figure out what > it is doing. Since I am using tri-linear element in 3D, I thought > the gradient of displacement could be computed with finite > difference (using center difference). However, it gives different > result than GetFEM++. So I am wondering how GetFEM++ computes the > gradients. > > best, > James > > > On Fri, May 11, 2012 at 5:24 AM, Andriy Andreykiv > <[email protected] <mailto:[email protected]>> > wrote: > > Dear James, > > When you speak about this linear ooperator, do you imply > element level or the global mesh level? > Basically, there are several ways to do this. > If you want to calculate strain in integration points (element > by element) please see calculation in > getfem_nonlinear_elasticty.h elasticity_nonlinear_term::compute. > > Below there is simple function, that extracts individual > linear strain components into a global vector that can be > postprocessed. > > The code: > > // > gefem::mesh m; // > getfem::mesh_fem mf_u(m); // the mesh_fem used to describe the > approximation of the displacements > set_classical_finite_element(approximation_order_U); > ... > > // compute displacement U from elastic problem > > .... > getfem::mesh_fem mf_strain(m); > > mf_strain.set_classical_discontinuous_finite_element(approximation_order_U > -1); > //approximation of strain has to be discontinous and one order > lower > plain_vector STRAIN_VECTOR(mf_strain.nb_dof()*6); /assuming 3D > with 6 components of strain in Voight notation > compute_strain(mf_u,mf_strain,U,STRAIN_VECTOR); > > // now STRAIN_VECTOR can be exported to VTK component by > component on mf_strain if you want to > > > ....... > //function compute_strain > void compute_strain(const getfem::mesh_fem& mf_u, const > getfem::mesh_fem& mf_strain, const plain_vector& U, > plain_vector& strain) > { > unsigned N = mf_u.linked_mesh().dim(); > GMM_ASSERT1(N==3,"This function is written as a demo only for > 3D"); > GMM_ASSERT1(gmm::vect_size(strain) == 6*mf_strain.nb_dof(), > "compute_strain: wrong size"); > > plain_vector GRAD(mf_strain.nb_dof() * N * N); > > //compute gradient of displacement on the whole mesh > getfem::compute_gradient(mf_u, mf_strain, U, GRAD); > > base_matrix E(N, N), gradU(N, N); > for (size_type i = 0; i < mf_strain.nb_dof(); i+=6) > { > //extract the gradient in the node > std::copy(GRAD.begin() + i * N*N, GRAD.begin()+(i + 1) * N*N, > gradU.begin()); > > // compute strain from the gradient > gmm::add(gradU,gmm::transposed(gradU),E); > gmm::scale(E, scalar_type(0.5)); > > //using Voight notation to convert the strain tensor > into a vector > // in 3D > strain[6*i ]=E(0,0); > strain[6*i+1]=E(1,1); > strain[6*i+2]=E(2,2); > strain[6*i+3]=E(0,1); > strain[6*i+4]=E(1,2); > strain[6*i+5]=E(0,2); > } > } > > > > > > > On 10 May 2012 21:26, James Zhou <[email protected] > <mailto:[email protected]>> wrote: > > Dear GetFEM++ users, > > In the context of solving a linear elasticity problem, is > there a way to extract the linear operator that converts > displacement into strain tensor? For example, if there > are n nodes, we would have a displacement vector of size > 3*n, if the strain tensors are reshaped into a 6*n vector > (assuming isotropic material), I would like to build the > 6*n x 3*n matrix that convert displacement into strain. > > I have tried a finite difference approach where strain is > written as the symmetric gradient of displacement, and the > gradient of displacement is approximated by finite > difference. However, this does not yield the same result > as the stress tensor computed in > src/getfem/getfem_derivatives.h > > Thanks! > > James > > _______________________________________________ > Getfem-users mailing list > [email protected] <mailto:[email protected]> > https://mail.gna.org/listinfo/getfem-users > > > > > > > > > _______________________________________________ > Getfem-users mailing list > [email protected] > https://mail.gna.org/listinfo/getfem-users -- Yves Renard ([email protected]) tel : (33) 04.72.43.87.08 Pole de Mathematiques, INSA-Lyon fax : (33) 04.72.43.85.29 20, rue Albert Einstein 69621 Villeurbanne Cedex, FRANCE http://math.univ-lyon1.fr/~renard ---------
_______________________________________________ Getfem-users mailing list [email protected] https://mail.gna.org/listinfo/getfem-users
