On 3/29/22 19:26, Léonhard YU wrote:
1) cell_mass_matrix(i, j) += rho *  fe_values.shape_value(i, q_point) *
   fe_values.shape_value(j, q_point) * fe_values.JxW(q_point);

Yes, this is not correct unless you prefix it with
  if (fe.system_to_component_index(i).first ==
      fe.system_to_component_index(j).first)


2) cell_mass_matrix(i, j) +=
         (mass_phi_i *  mass_phi_j) *
          fe_values.JxW(q_point) * rho;

     where  mass_phi_i:
             const Tensor<2, dim>
                      mass_phi_i  =get_mass_matrix(fe_values, i, q_point);

             (using this function)
   template <int dim>
   inline Tensor<2, dim>  get_mass_matrix(const FEValues<dim> &fe_values,
       const unsigned int shape_func_i,
       const unsigned int shape_func_j,
       const unsigned int q_point)
{
   Tensor<2, dim> tmp;
   for (unsigned int t = 0; t < dim; ++t)
       for (unsigned int k = 0; k < dim; ++k)
       {
          tmp[t][k] = fe_values.shape_value_component(shape_func_i, q_point, t)*                   fe_values.shape_value_component(shape_func_j, q_point, k);
       }
   return tmp;
}

I don't understand this code. You say
  const Tensor<2, dim> mass_phi_i
       = get_mass_matrix(fe_values, i, q_point);

but the definition of get_mass_matrix() takes four arguments.


2) cell_mass_matrix(i, j) +=
         (mass_phi_i *  mass_phi_j) *
          fe_values.JxW(q_point) * rho;

     where  mass_phi_i:
             const  SymmetricTensor  <2, dim>
                      mass_phi_i  =get_mass_matrix(fe_values, i, q_point);

             (using this function)
   template <int dim>
   inline SymmetricTensor<2, dim>
    get_mass_matrix  (const FEValues<dim> &fe_values,
              const unsigned int   shape_func,
              const unsigned int   q_point)
   {
     SymmetricTensor<2, dim> tmp;
     for (unsigned int i = 0; i < dim; ++i)
       tmp[i][i] = fe_values.shape_value_component(shape_func, q_point, i) *
                        fe_values.shape_value_component(shape_func, q_point, i);

     for (unsigned int i = 0; i < dim; ++i)
       for (unsigned int j = i + 1; j < dim; ++j)
         tmp[i][j] =
           fe_values.shape_value_component(shape_func, q_point, i) *
            fe_values.shape_value_component(shape_func, q_point, j);
     return tmp;
   }

This suffers from the same problem.

I think you want to read through the module on vector-valued problems to understand the general philosophy about dealing with vector problems.

Best
 W.

--
------------------------------------------------------------------------
Wolfgang Bangerth          email:                 [email protected]
                           www: http://www.math.colostate.edu/~bangerth/

--
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].
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/42ff523d-ecae-e8dd-6417-2cf620301c87%40colostate.edu.

Reply via email to