Hi Wolfgang and Jörg I am taking advantage for a second question. I have made a few minor changes in the code and re-implemented it. In serial version, all works fine so far. However, when running in parallel, I am seeing an issue in the method PlasticityContactProblem::update_solution_and_constraints.
In particular, it turns out that the value of const unsigned int index_z = dof_indices[q_point]; might be out of the range of TrilinosWrappers::MPI::Vector <http://www.dealii.org/8.5.0/doxygen/deal.II/classTrilinosWrappers_1_1MPI_1_1Vector.html> lambda(locally_relevant_dofs, mpi_communicator); lambda = newton_rhs_uncondensed; I have solved run-time issues by checking if if ( lambda.in_local_range(index_z)) but I wonder if the issue is possible (in principle, I mean) for the published code or is it my implementation that is not correct? Note that I am basically implementing step-42 upon the code built in step-20 using PETScWrappers::MPI::Vector, rather than step-40. If of use, I am writing my method below. Thanks for your help. Alberto ::update_solution_and_constraints () // The following function is the first function we call in each outer loop of the Newton iteration. // What it does is to project the solution onto the feasible set and update the active set for the degrees of freedom that touch or penetrate the obstacle. { std::pair <unsigned, unsigned> vrange; this->IO.log() << " In update_solution_and_constraints. " << std::flush; // We need to write into the solution vector (which we can only do with fully distributed vectors without ghost elements) and we need to read the Lagrange multiplier and the elements of the diagonal mass matrix from their respective vectors (which we can only do with vectors that do have ghost elements), so we create the respective vectors. We then also initialize the constraints object that will contain constraints from contact, as well as an object that contains an index set of all locally owned degrees of freedom that are part of the contact: std::vector<bool> dof_touched( this->dof_handler.n_dofs(), false); PETScWrappers::MPI::Vector distributed_solution( this->locally_owned_dofs, this->mpi_communicator ); distributed_solution = this->accumulated_displacement; vrange = distributed_solution.local_range(); this->IO.log() << " distributed_solution local range = " << vrange.first << ", " << vrange.second << "\n" << std::flush; PETScWrappers::MPI::Vector lambda( this->locally_relevant_dofs, this ->mpi_communicator); lambda = this->system_rhs_uncondensed; vrange = lambda.local_range(); this->IO.log() << " lambda local range = " << vrange.first << ", " << vrange.second << "\n" << std::flush; PETScWrappers::MPI::Vector diag_mass_matrix_vector_relevant( this->locally_relevant_dofs, this->mpi_communicator); diag_mass_matrix_vector_relevant = this->diag_mass_matrix_vector; vrange = diag_mass_matrix_vector_relevant.local_range(); this->IO.log() << " diag_mass_matrix_vector_relevant local range = " << vrange.first << ", " << vrange.second << "\n" << std::flush; this->all_constraints.reinit( this->locally_relevant_dofs ); this->contact_constraints.reinit( this->locally_relevant_dofs ); this->active_set.clear(); // The second part is a loop over all cells in which we look at each point where a degree of freedom is defined whether the active set condition is true and we need to add this degree of freedom to the active set of contact nodes. As we always do, if we want to evaluate functions at individual points, we do this with an FEValues object (or, here, an FEFaceValues object since we need to check contact at the surface) with an appropriately chosen quadrature object. We create this face quadrature object by choosing the "support points" of the shape functions defined on the faces of cells (for more on support points, see this glossary entry). As a consequence, we have as many quadrature points as there are shape functions per face and looping over quadrature points is equivalent to looping over shape functions defined on a face. With this, the code looks as follows: Quadrature<dim-1> face_quadrature( this ->fe.get_unit_face_support_points()); FEFaceValues<dim> fe_values_face( this->fe, face_quadrature, update_quadrature_points ); const unsigned int dofs_per_face = this->fe.dofs_per_face; const unsigned int n_face_q_points = face_quadrature.size(); std::vector<types::global_dof_index> dof_indices(dofs_per_face); typename DoFHandler<dim>::active_cell_iterator cell = this->dof_handler.begin_active(), endc = this->dof_handler.end(); for (; cell != endc; ++cell) if (!cell->is_artificial()) { if (cell->is_locally_owned()) // MY OWN CODE { for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face) // the Contact boundary is set to value 3, i.e. set_boundary_id (3); if (cell->face(face)->at_boundary() && cell->face(face)->boundary_id() == 3) { fe_values_face.reinit(cell, face); cell->face(face)->get_dof_indices(dof_indices); for (unsigned int q_point=0; q_point<n_face_q_points; ++q_point) { // At each quadrature point (i.e., at each support point of a degree of freedom located on the contact boundary), we then ask whether it is part of the X2-displacement degrees of freedom and if we haven't encountered this degree of freedom yet (which can happen for those on the edges between faces), we need to evaluate the gap between the deformed object and the obstacle. If the active set condition is true, then we add a constraint to the ConstraintMatrix object that the next Newton update needs to satisfy, set the solution vector's corresponding element to the correct value, and add the index to the IndexSet object that stores which degree of freedom is part of the contact: const unsigned int component = this ->fe.face_system_to_component_index(q_point).first; const unsigned int index_X2 = dof_indices[q_point]; // MY OWN CODE: component has been changed to 1 since the direction of contact is X2 for the cell. if ((component == 1) && (dof_touched[index_X2] == false)) { dof_touched[index_X2] = true; const Point<dim> this_support_point = fe_values_face.quadrature_point(q_point); const double obstacle_value = this->obstacle->value(this_support_point, 1); const double solution_here = this->accumulated_displacement(index_X2); // solution(index_X2); // I believe it should be undeformed_gap >=0, i.e. : const double undeformed_gap = this_support_point(1) - obstacle_value; // rather than obstacle_value - this_support_point(1); double e_modulus; if ( cell->material_id() == 0 ) e_modulus = this ->cytoplasm_model.youngmodulus(); else if ( cell->material_id() == 1 ) e_modulus = this ->nucleus_model.youngmodulus(); else e_modulus = 1; if ( lambda.in_local_range(index_X2) && diag_mass_matrix_vector_relevant.in_local_range(index_X2) ) { lambda(index_X2); diag_mass_matrix_vector_relevant(index_X2); this->hanging_node_constraints.is_constrained(index_X2); const double c = 100.0 * e_modulus; if ((lambda(index_X2) / diag_mass_matrix_vector_relevant(index_X2) + c * ( - solution_here - undeformed_gap ) // I believe it should be this, rather than (solution_here - undeformed_gap) > 0) && !(this->hanging_node_constraints.is_constrained(index_X2)) ) { this->contact_constraints.add_line( index_X2 ); // it was this->all_constraints.add_line(index_X2); this->contact_constraints.set_inhomogeneity(index_X2, - undeformed_gap); // it was this->all_constraints.set_inhomogeneity // I believe it should be this, rather than undeformed_gap); distributed_solution( index_X2 ) = - undeformed_gap; // I believe it should be this, rather than undeformed_gap; this->active_set.add_index( index_X2 ); } }; } } } } } // At the end of this function, we exchange data between processors updating those ghost elements in the solution variable that have been written by other processors. distributed_solution.compress(VectorOperation::insert); this->accumulated_displacement = distributed_solution; // We then merge the constraints from hanging nodes into the ConstraintMatrix object that already contains the active set. this->contact_constraints.close(); // it was this->all_constraints.close(); this->all_constraints.merge(this->hanging_node_constraints); this->all_constraints.merge(this->contact_constraints); this->pcout << " Size of active set: " << Utilities::MPI::sum((this->active_set & this->locally_owned_dofs).n_elements(), this->mpi_communicator) << std::endl; } *Alberto Salvadori* Dipartimento di Ingegneria Civile, Architettura, Territorio, Ambiente e di Matematica (DICATAM) Universita` di Brescia, via Branze 43, 25123 Brescia Italy tel 030 3711239 fax 030 3711312 e-mail: alberto.salvad...@unibs.it web-pages: http://m4lab.unibs.it/faculty.html http://dicata.ing.unibs.it/salvadori On Wed, Jan 17, 2018 at 11:01 PM, Frohne, Joerg < fro...@mathematik.uni-siegen.de> wrote: > >>On 01/17/2018 02:45 PM, Frohne, Joerg wrote: > >> obviously this is a mistake. I checked the source files which we have > used for results in corresponding paper. > >> There I have found the following coding: > >> > >> resid_old = resid; > >> > >> resid_vector = system_rhs_newton; > >> resid_vector.compress(VectorOperation::insert); > >> > >> int is_my_set_changed = (active_set == active_set_old) ? 0 : > 1; > >> int num_changed = Utilities::MPI::sum(is_my_set_changed, > >> MPI_COMM_WORLD); > >> if (num_changed == 0 && resid < 1e-8) > >> break; > >> active_set_old = active_set; > >> > >> Unfortunately I have not been working with deal.ii for over two years > now. > >> For this reason I would be glad if someone else could fix this bug. > >> Otherwise I would need a short manual how to do it by myself. > > >I can do that for you. Do I understand correctly that the only change > >that needs to happen is to move > > active_set_old = active_set; > >to *after* the if-statement? > > Exactly. > > Best, > Joerg > > -- > 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. > -- Informativa sulla Privacy: http://www.unibs.it/node/8155 -- 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.
