Indeed, one commonly for stress recovery used approach is nodal averaging. See for example Zienkiewicz and Zhu "A Simple Estimator and Adaptive Procedure for Practical Engineering Analysis" http://onlinelibrary.wiley.com/doi/10.1002/nme.1620240206/abstract You can obtain the best fit in L2 norm by doing an approximate inversion of the mass matrix, i.e. using lumped mass matrix. I think it is a good compromise between a proper L2 projection which might be slow for big meshes and heuristics of averaging over neighbouring elements.
Regards, Denis. On Saturday, October 31, 2015 at 8:34:59 PM UTC+1, Timo Heister wrote: > > Nick, > > I can not think of a much simpler way of doing what you are trying to > do and using the global dof index is a reasonable way to do this. > > Two comments: > - Your code is dimension dependent, but I am sure you are aware of this. > :-) > - The stress/strain is discontinuous on the interface between cells so > your problem is not well-defined and taking an average like you do, is > not the "right" approach in the FEM world. An alternative would be to > define a FiniteElement (for example a continuous element) and then > project the stress into that space (basically finding the > best-approximation). This involves assembling and solving a global > mass matrix. Armed with this field, you can evaluate it where ever you > want. > > Best, > Timo > > On Sat, Oct 31, 2015 at 9:45 AM, Nicholas Padon <[email protected] > <javascript:>> wrote: > > Hello, > > > > I'm solving a 2-D plane strain elasticity problem in line with step-8 > (and > > some of the techniques from step-18), and I'm trying to output the > stresses > > both at the cell centers and at the support points (for arbitrary > Lagrange > > elements). I realize that the stresses at the support points are > > discontinuous between cells, so I've also got to average them in some > way. > > I've gotten it working (I think), but I'm wondering if there's a better > way > > to do what I want. Basically I've: > > > > 1. Created a quadrature formula based on the support points of a > particular > > cell > > 2. On each cell, loop through the support points and add up some total > > strain, which is later multiplied by C_ijkl and divided by the number of > > support points (to get the average stress per cell) > > 3. On each cell, add each stress component to a vector that keeps track > of > > support point dofs > > 4. Outside the cell loop, divide each support point stress component by > the > > number of times that the support point was written to (to get the > "average" > > stress at the support point) > > > > In particular, I'm wondering if there's a more clever way to store the > dofs > > associated with each support point. Right now I'm writing the stress > > component to the 0th dof at each support point because I couldn't figure > out > > a better way to associate the entries in, say, sigma_xx_nodes, with the > node > > ordering that DataOut expects. Any thoughts on improving the method are > > appreciated... > > > > My code (imagine in an ::output_results() function) is: > > > > const std::vector<Point<dim>> support_points = > > fe.base_element(0).get_unit_support_points(); > > std::vector<double> weights(support_points.size(),1); > > const Quadrature<dim> node_quadrature_formula(support_points,weights); > > > > double num_cells = triangulation.n_active_cells(); > > types::global_dof_index n_dofs = dof_handler.n_dofs(); > > unsigned int n_node_q_points = node_quadrature_formula.size(); > > double num_vertices = triangulation.n_used_vertices(); > > > > std::cout << "Num cells: " << num_cells << std::endl; > > std::cout << "Num vertices: " << num_vertices << std::endl; > > std::cout << "Total DOF: " << n_dofs << std::endl; > > std::cout << "Length of solution vector:" << solution.size() << > std::endl; > > std::cout << "Number of support points: " << support_points.size() << > > std::endl; > > std::cout << "Number of nodal quad points:" << > > node_quadrature_formula.size() << std::endl; > > > > > > FEValues<dim> fe_node_values (fe, node_quadrature_formula, > > update_values | update_gradients | > > update_quadrature_points | > update_JxW_values); > > > > typename DoFHandler<dim>::active_cell_iterator cell = > > dof_handler.begin_active(), > > endc = > dof_handler.end(); > > > > std::vector< SymmetricTensor<2,dim>> > node_strain_results(n_node_q_points); > > > > SymmetricTensor<2,dim> cell_strain; > > SymmetricTensor<2,dim> cell_stress; > > > > Vector<double> sigma_xx_cell(num_cells); > > Vector<double> sigma_yy_cell(num_cells); > > Vector<double> sigma_xy_cell(num_cells); > > > > std::vector<int> number_of_stresses (n_dofs,0); > > Vector<double> sigma_xx_nodes(n_dofs); > > Vector<double> sigma_yy_nodes(n_dofs); > > Vector<double> sigma_xy_nodes(n_dofs); > > > > std::vector<types::global_dof_index> > global_dof_indices(fe.dofs_per_cell); > > > > int cell_index =0; > > const FEValuesExtractors::Vector displacements (0); > > for(;cell!=endc;++cell) > > { > > fe_node_values.reinit(cell); > > cell_strain.clear(); > > > > //get the strains at each support point > > > fe_node_values[displacements].get_function_symmetric_gradients(solution,node_strain_results); > > > > > > //compute the average stresses at each cell first > > cell->get_dof_indices(global_dof_indices); > > for(int q=0;q<n_node_q_points;++q) > > { > > cell_stress.clear(); > > //compute the stress at the support point > > cell_stress = (stress_strain_tensor)*node_strain_results[q]; > > > > //add the stress at the support point to the correct dof index > > //NOTE: because support points can be shared between cells, these > > //vectors will need to be averaged based on the number of times they've > > //been written to > > sigma_xx_nodes(global_dof_indices[q*dim]) += cell_stress[0][0]; > > sigma_yy_nodes(global_dof_indices[q*dim]) += cell_stress[1][1]; > > sigma_xy_nodes(global_dof_indices[q*dim]) += cell_stress[0][1]; > > number_of_stresses[global_dof_indices[q*dim]] += 1; > > > > //add the strains up at each evaluation point > > cell_strain += node_strain_results[q]; > > } > > > > //multiply the stress/strain tensor divide by the number of quadrature > > points > > //to get the average stress in the cell > > cell_stress.clear(); > > cell_stress = (stress_strain_tensor*cell_strain)/n_node_q_points; > > //assign each stress to the right cell > > sigma_xx_cell[cell_index] = cell_stress[0][0]; > > sigma_yy_cell[cell_index] = cell_stress[1][1]; > > sigma_xy_cell[cell_index] = cell_stress[0][1]; > > cell_index++; > > } > > for(types::global_dof_index i=0;i<n_dofs;++i) > > { > > if(number_of_stresses[i]>0) > > { > > sigma_xx_nodes(i) = sigma_xx_nodes(i)/number_of_stresses[i]; > > sigma_yy_nodes(i) = sigma_yy_nodes(i)/number_of_stresses[i]; > > sigma_xy_nodes(i) = sigma_xy_nodes(i)/number_of_stresses[i]; > > } > > } > > > > data_out.add_data_vector(sigma_xx_cell,"sigma_xx_cell_center"); > > data_out.add_data_vector(sigma_yy_cell,"sigma_yy_cell_center"); > > data_out.add_data_vector(sigma_xy_cell,"sigma_xy_cell_center"); > > data_out.add_data_vector(sigma_xx_nodes,"sigma_xx_nodes"); > > data_out.add_data_vector(sigma_yy_nodes,"sigma_yy_nodes"); > > data_out.add_data_vector(sigma_xy_nodes,"sigma_xy_nodes"); > > > > > > Thanks in advance. > > > > Nick Padon > > > > -- > > 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 [email protected] <javascript:>. > > For more options, visit https://groups.google.com/d/optout. > > > > -- > Timo Heister > http://www.math.clemson.edu/~heister/ > -- 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 [email protected]. For more options, visit https://groups.google.com/d/optout.
