Dear Alberto, In my opinion both of the approaches that you've outlined are plausible ways of implementing Dirichlet constraints, but there are a couple of key differences that I can quickly outline. I'll refer to the approach listed in your first post as method 1 and the split approach in your second post as method 2.
Method 1: - You need to create a union of component masks for the field components to be constrained. You can do this using ComponentMask::operator| <https://www.dealii.org/8.5.0/doxygen/deal.II/classComponentMask.html#a9070d2c18038c413113775e420cf9fef>, and there are examples of its use in the tutorials. So presumably > ComponentMask displacements_mask = fe.component_mask (displacements); would change to something like > const ComponentMask displacements_and_concentration_mask = > fe.component_mask (displacements) | fe.component_mask (concentration); - You need to implement the definition of all constraints for a boundary into a single class, i.e. for the function > IncrementalDirichletBoundaryValues<dim>::vector_value (const Point<dim> & > p , Vector<double> & return_value) const all of the constrained components (as selected via the component mask) of return_value need to be sensibly defined. - With this monolithic approach (and unless you do something strange), you should never have conflicting definitions of a constraint on a boundary. Method 2: - You are first imposing one set of constraints, and then a second. This is now a sequential operation rather than a monolithic one as used in method 1. - You need only have one std::map<types::global_dof_index,double> boundary_values, sequentially interpolate all of the Dirichlet constraint definitions into the map and apply them once. This is more computationally efficient (you only modify the linear system once). - Applying multiple boundary conditions to subsets of DoFs on a boundary introduces the possibility of overwriting constraints (e.g. by accidentally selecting the same component of a ComponentMask twice). So you could introduce a bug where the second constraint definition dominates the first. There is a warning on this point in the documentation to VectorTools::interpolate_boundary_values <https://www.dealii.org/8.5.0/doxygen/deal.II/namespaceVectorTools.html#a9f3e3ae1396811f998cc35f94cbaa926>. This would normally only require consideration for DoFs shared between adjacent Dirichlet boundaries, but now some extra care must be taken in this situation. If it makes any difference, my personal preference is to use method 2. With this approach you can define interesting boundary conditions once and then easily mix and match which are applied to which components of different boundaries. But, of course, opinions would differ based on design philosophies and (perhaps) mathematical rigour related to the theory on boundary value problems. Best, J-P On Thursday, June 1, 2017 at 1:54:23 AM UTC+2, Alberto Salvadori wrote: > > I am adding here an attempt I made. It seems to work but since this was > more intuition rather than full understanding, I do appreciate your > comments. > So, this is what I did: basically, I created two Dirichlet boundary > conditions and two masks, and I applied conditions in sequence, like this: > > > // Dirichlet bcs > > // ------------- > > > > PETScWrappers::MPI::Vector tmp (locally_owned_dofs,mpi_communicator); > > > // Dirichlet bc for displacements > > { > > std::map<types::global_dof_index,double> boundary_values; > > > ComponentMask displacements_mask = fe.component_mask (displacements); > > VectorTools::interpolate_boundary_values (dof_handler, > > 0, > > > IncrementalDirichletBoundaryValues<dim>( TIDM, NRit, GPDofs ), > > boundary_values, > > displacements_mask); > > > > MatrixTools::apply_boundary_values (boundary_values, > > system_matrix, > > tmp, > > system_rhs, > > false); > > } > > > // Dirichlet bc for concentrations > > { > > std::map<types::global_dof_index,double> boundary_values; > > > > ComponentMask displacements_mask = fe.component_mask (concentration); > > VectorTools::interpolate_boundary_values (dof_handler, > > 0, > > > IncrementalDirichletBoundaryValuesForConcentration<dim>( TIDM, NRit, > GPDofs ), > > boundary_values, > > displacements_mask); > > > > PETScWrappers::MPI::Vector tmp (locally_owned_dofs,mpi_communicator); > > MatrixTools::apply_boundary_values (boundary_values, > > system_matrix, > > tmp, > > system_rhs, > > false); > > } > > > incremental_displacement = tmp; > > > Within the classes IncrementalDirichletBoundaryValuesForConcentration > I have returned something like: > > template <int dim> > > void > > IncrementalDirichletBoundaryValuesForConcentration<dim>:: > > vector_value (const Point<dim> & p , > > Vector<double> & return_value) const > > { > > > > static const unsigned int c_component = dim + 2; > > return_value( c_component ) = ... ; > > } > > > Does it make sense? > > Thanks > Alberto > > Il giorno mercoledì 31 maggio 2017 17:48:13 UTC-4, Alberto Salvadori ha > scritto: >> >> Dear all, >> >> your help is appreciated about how component mask and Dirichlet bc >> application work. >> >> I am implementing a SmallStrainDiffusionMechanicalProblem class, with 4 >> fields: displacements, pressure, dilatation, and concentration. >> >> template <int dim> >> >> class SmallStrainDiffusionMechanicalProblem >> >> { >> >> ... >> >> >> // dofs definiton and dofs block enumeration >> >> // dim = displacements dofs >> >> // 1 = p >> >> // 1 = J >> >> // 1 = c ( interstitial concentration ) >> >> const unsigned int GPDofs = dim + 3; >> >> >> static const unsigned int first_u_component = 0; >> >> static const unsigned int p_component = dim; >> >> static const unsigned int J_component = dim + 1; >> >> static const unsigned int c_component = dim + 2; >> >> >> enum >> >> { >> >> u_dof = 0, // displacement block ( dim components ) >> >> p_dof = 1, // pressure block ( one component ) >> >> J_dof = 2, // dilatation block ( one component ) >> >> c_dof = 3 // concentration block ( one component ) >> >> }; >> >> >> >> ... >> >> } >> >> >> In a former implementation that did not include concentration fields, in >> order to impose bc to the vector solution incremental_displacement, I >> have edited some code from the tutorials, like this: >> >> >> const FEValuesExtractors::Vector displacements (first_u_component); >> >> const FEValuesExtractors::Scalar pressure(p_component); >> >> const FEValuesExtractors::Scalar dilatation(J_component); >> >> const FEValuesExtractors::Scalar concentration(c_component); >> >> ...... >> >> std::map<types::global_dof_index,double> boundary_values; >> >> >> ComponentMask displacements_mask = fe.component_mask (displacements); >> >> VectorTools::interpolate_boundary_values (dof_handler, >> >> 0, >> >> >> IncrementalDirichletBoundaryValues<dim>( TIDM, NRit, GPDofs ), >> >> boundary_values, >> >> displacements_mask); >> >> >> >> PETScWrappers::MPI::Vector tmp (locally_owned_dofs,mpi_communicator); >> >> MatrixTools::apply_boundary_values (boundary_values, >> >> system_matrix, >> >> tmp, >> >> system_rhs, >> >> false); >> >> incremental_displacement = tmp; >> >> I would now impose Dirichlet bc on both displacements and concentrations >> and assume that the Dirichlet boundary for both fields is the same. Shall I >> therefore change the mask, I assume. Can it be done like this? >> >> ComponentMask displacements_mask = fe.component_mask (displacements, >> concentrations); >> >> or shall two masks be defined or even else? >> >> Many thanks! >> Alberto >> >> >> -- 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.
