Dear all, I am trying to port the previous linear elastic code to MeshWorker, but I have some questions about this class.
First of all, I have read this question <https://groups.google.com/forum/#!searchin/dealii/meshworker|sort:relevance/dealii/JTIM1zASACU/wIvwrU2hAQAJ>, and I read that MeshWorker is sort of abandoned: the issue with MeshWorker is that when it was implemented, we didn't write > enough documentation/tutorial programs/examples [...]. Development then > also stopped, and so [...] nobody today knows any more whether they can or > can not be done with MeshWorker. The only person who still knows how this > all works unfortunately > also doesn't regularly participate in the help forum any more :-( > > [...] you can of course ask questions about how this or that can be done, > but you're likely not going to get a lot of helpful answers because nobody > knows any more how this works. It's really an unfortunate situation we got > ourselves into with this small corner of the library. > So my first question is, should I avoid using this class and implement parallel loops by hand (via TBB or other means)? Then, my doubts about the use of MeshWorkers. I cannot really grasp it even after reading tutorials 12, 16, and 39. What is the purpose of the cell, boundary, and face functions? The cell is called exactly on each n-dimensional cell, I've checked that number, what about the other two? Boundary is called, as I read, for each n-dimensional cell on the boundary (numbers again match). Face function is called on each internal (n-1)-dimensional face (not on the boundary), as far as I understand. Is this correct? When assigning a Neumann condition on faces, for instance in elasticity apply a load, where is the correct place? I have tried to get the code into the cell function, but I didn't find a way to get the faces, boundary markers, and the face values. My code is at the bottom, if you need it. Is there an analogous of the FullMatrix for the right hand side vector? The only thing similar to a vector in the DoFInfo is a BlockVector... Thanks for any clarification you can give me. Cheers! // Integrate local cell matrix void integrate_cell_term(MeshWorker::DoFInfo<dim> &dinfo, MeshWorker:: IntegrationInfo<dim> &info) { loopc++; // FE Values const FEValuesBase<dim> &fe_v = info.fe_values(); // Local Lame values, here are stupid constants std::vector<double> lambda_values(fe_v.n_quadrature_points), mu_values(fe_v.n_quadrature_points); // Local element matrix FullMatrix<double> &local_matrix = dinfo.matrix(0).matrix; // Local element RHS vector BlockVector<double> &local_vector = dinfo.vector(0); // Jacobian const std::vector<double> &JxW = fe_v.get_JxW_values (); // Local Lame values lambda.value_list(fe_v.get_quadrature_points(), lambda_values); mu.value_list(fe_v.get_quadrature_points(), mu_values); // For each DOF of a cell (i-th) for (unsigned int i = 0; i < fe_v.dofs_per_cell; ++i) { const unsigned int component_i = info.fe_values().get_fe(). system_to_component_index(i).first; // For each DOF of a cell (j-th) for (unsigned int j = 0; j < fe_v.dofs_per_cell; ++j) { const unsigned int component_j = info.fe_values().get_fe(). system_to_component_index(j).first; // Compute the (i, j)-th component for LHS via quadrature for (unsigned int q_point = 0; q_point < fe_v. n_quadrature_points; ++q_point) { local_matrix(i,j) += ((fe_v.shape_grad(i, q_point)[ component_i] * fe_v.shape_grad(j, q_point)[ component_j] * lambda_values[q_point]) + (fe_v.shape_grad(i, q_point)[ component_j] * fe_v.shape_grad(j, q_point)[ component_i] * mu_values[q_point]) + ( (i == j) ? (fe_v.shape_grad(component_i, q_point) * fe_v.shape_grad(component_j, q_point) * mu_values[q_point]) : 0) ) * JxW[q_point]; } // Integral for bilinear form } // j // Computed values for the RHS std::vector<Vector<double>> rhs_values(fe_v.n_quadrature_points, Vector<double>(dim)); // Computed values for the RHS std::vector<Vector<double>> rhs_face_values(fe_v.n_quadrature_points, Vector<double>(dim)); // RHS values right_hand_side.vector_value_list(fe_v.get_quadrature_points(), rhs_values); // Compute the i-th component for RHS via quadrature (body forces) for(unsigned int q_point = 0; q_point < fe_v.n_quadrature_points; ++q_point) { local_vector(i) += fe_v.shape_value(i, q_point) * rhs_values[q_point](component_i) * JxW[q_point]; } // Integral for linear form } // i // Computed values for the RHS on the faces std::vector<Vector<double>> rhs_face_values(fe_v.n_quadrature_points, Vector<double>(dim)); // Neumann (boundary forces on top boundary == id 2) for (unsigned int face_number = 0; face_number < GeometryInfo<dim>:: faces_per_cell; ++face_number) { loopf++; if (dinfo.cell->face(face_number)->at_boundary() && (dinfo.cell-> face(face_number)->boundary_id() == 2)) { // FE Values on face (???) const FEValuesBase<dim> &fe_face_values = info.fe_values(); // Get face FE values for this face and cell fe_face_values.reinit(dinfo.cell, face_number); // Neumann force vector on the face right_face_side.vector_value_list(fe_face_values. get_quadrature_points(), rhs_face_values); for(unsigned int f_q_point = 0; f_q_point < fe_v. n_quadrature_points; ++f_q_point) { // Add contributions to all DOFs for (unsigned int i = 0; i < fe_v.dofs_per_cell; ++i) { const unsigned int component_i = fe_v.system_to_component_index(i).first; local_vector(i) += (fe_face_values.shape_value(i, f_q_point) * rhs_face_values[f_q_point][component_i]) * fe_face_values.JxW(f_q_point); } } // Integral for linear form } } // Boundary forces } -- 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.
