Hi, I found some issue when constructing the FE B-operator for multiple elements. In a single element it works perfectly. My approach is the following:
template<int dim> FullMatrix<double> get_b_operator (const FEValues<dim> &fe_values, const unsigned int dofs_per_cell, const unsigned int q) { FullMatrix<double> tmp(dim, GeometryInfo<dim>::vertices_per_cell); std::vector<DerivativeForm<1, dim, dim> > invJ = fe_values.get_inverse_jacobians(); for (unsigned int i = 0; i < dofs_per_cell; i = i + dim) { const unsigned int index = i / dim; // COMPUTE: B-operator (Remark: This version has to be extended for 3D!) tmp[0][index] = invJ[q][0][0] * fe_values.shape_grad(i, q)[0] + invJ[q][0][1] * fe_values.shape_grad(i, q)[1]; tmp[1][index] = invJ[q][1][0] * fe_values.shape_grad(i, q)[0] + invJ[q][1][1] * fe_values.shape_grad(i, q)[1]; } return tmp; } Hence, I am computing B = inv(J) * transpose(dN). For a single element I receive the following correct results for dN: *C++ CODE* SHAPE FUNCTION DERIVATIVES (GAUSS POINT 1) -0.788675 -0.788675 0.788675 -0.211325 0.211325 0.211325 -0.211325 0.788675 SHAPE FUNCTION DERIVATIVES (GAUSS POINT 2) -0.788675 -0.211325 0.788675 -0.788675 0.211325 0.788675 -0.211325 0.211325 SHAPE FUNCTION DERIVATIVES (GAUSS POINT 3) -0.211325 -0.788675 0.211325 -0.211325 0.788675 0.211325 -0.788675 0.788675 SHAPE FUNCTION DERIVATIVES (GAUSS POINT 4) -0.211325 -0.211325 0.211325 -0.788675 0.788675 0.788675 -0.788675 0.211325 *DEAL.II* SHAPE FUNCTION DERIVATIVES -0.788675 -0.788675 -0.788675 -0.788675 0.788675 -0.211325 0.788675 -0.211325 -0.211325 0.788675 -0.211325 0.788675 0.211325 0.211325 0.211325 0.211325 -0.788675 -0.211325 -0.788675 -0.211325 0.788675 -0.788675 0.788675 -0.788675 -0.211325 0.211325 -0.211325 0.211325 0.211325 0.788675 0.211325 0.788675 -0.211325 -0.788675 -0.211325 -0.788675 0.211325 -0.211325 0.211325 -0.211325 -0.788675 0.788675 -0.788675 0.788675 0.788675 0.211325 0.788675 0.211325 -0.211325 -0.211325 -0.211325 -0.211325 0.211325 -0.788675 0.211325 -0.788675 -0.788675 0.211325 -0.788675 0.211325 0.788675 0.788675 0.788675 0.788675 Due to the ordering of the element connectivity (numbering of nodes are different in deal.II) we have to swap column 3 and 4 for comparison. Now for 4 elements there is something strange when I compute my shape function derivatives within deal.II: *DEAL.II* SHAPE FUNCTION DERIVATIVES -1.577350 -1.577350 -1.577350 -1.577350 1.577350 -0.422650 1.577350 -0.422650 -0.422650 1.577350 -0.422650 1.577350 0.422650 0.422650 0.422650 0.422650 -1.577350 -0.422650 -1.577350 -0.422650 1.577350 -1.577350 1.577350 -1.577350 -0.422650 0.422650 -0.422650 0.422650 0.422650 1.577350 0.422650 1.577350 -0.422650 -1.577350 -0.422650 -1.577350 0.422650 -0.422650 0.422650 -0.422650 -1.577350 1.577350 -1.577350 1.577350 1.577350 0.422650 1.577350 0.422650 -0.422650 -0.422650 -0.422650 -0.422650 0.422650 -1.577350 0.422650 -1.577350 -1.577350 0.422650 -1.577350 0.422650 1.577350 1.577350 1.577350 1.577350 The results for one element change now by a factor of 2 and with increasing element size it always differs by another factor of 2 again. I was wondering, if my approach to compute my B-operator was wrong. Doesn't fe_values.shape_grad() give me the shape function derivatives directly? Or do I have to use shape_grad_component? Additionally, I checked the symmetric_gradient from step-18 and it contains actually the correct values I need, but somehow in a strange order: SYMMETRIC GRADIENT (GAUSS POINT 1) 0.000000 0.211325 0.211325 0.422650 SYMMETRIC GRADIENT (GAUSS POINT 2) 0.000000 0.211325 0.211325 1.577350 SYMMETRIC GRADIENT (GAUSS POINT 3) 0.000000 0.788675 0.788675 0.422650 SYMMETRIC GRADIENT (GAUSS POINT 4) 0.000000 0.788675 0.788675 1.577350 Kind regards, S. A. Mohseni -- The deal.II project is located at http://www.dealii.org/ For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en --- You received this message because you are subscribed to the Google Groups "deal.II User Group" group. To unsubscribe from this group and stop receiving emails from it, send an email to dealii+unsubscr...@googlegroups.com. For more options, visit https://groups.google.com/d/optout.