Dirichlet bc via apply_boundary_values function

[Мария Бронзова]

Dear all,

I"ve been trying to impose pressure boundary conditions in the context of a mixed displacement-pressure formulation problem, applying the "apply_boundary_values" function. I am dealing with complex numbers in my case, so I am imposing complex pressure values on certain faces with boundary id 1, as defined further:

 template <int dim>

  class PressureBoundaryValues : public Function<dim, std::complex<double>>

  {

  public:

  PressureBoundaryValues()

    : Function<dim, std::complex<double>>(dim+1)

      {}

  virtual std::complex<double> value(const Point<dim> &p,

         const unsigned int component = 0) const override;

  };

  template <int dim>

  std::complex<double> PressureBoundaryValues<dim>::value(const Point<dim> & p,

               const unsigned int component) const

  {

  const std::complex<double> i = {0,1};

  return 250.+ 0.0*i;

  }

for (const auto &cell : triangulation.active_cell_iterators())

for (const auto &face : cell->face_iterators())

if ((face->center()[dim-1] == 0))

face->set_all_boundary_ids(1);

FEValuesExtractors::Scalar pressure(dim);

FEValuesExtractors::Vector velocities(0);

std::map<types::global_dof_index, std::complex<double>> boundary_values;

VectorTools::interpolate_boundary_values(dof_handler,

    1,

     PressureBoundaryValues<dim>(),

     boundary_values,

     fe.component_mask(pressure));

MatrixTools::apply_boundary_values(boundary_values,

  system_matrix,

  solution,

  system_rhs);

The problem is that the real part gets perfectly assigned, whereas the imaginary value is ignored completely, even though I was trying to prescribe a zero value. In the results of calculation i get everything possible, rather than zero.

Could someone assist me in this case, please?

Thank you a lot!

With kind regards,

Mariia

[Wolfgang Bangerth]

On 10/12/21 7:31 AM, Мария Бронзова wrote:
>
> The problem is that the real part gets perfectly assigned, whereas the
> imaginary value is ignored completely, even though I was trying to
> prescribe a zero value. In the results of calculation i get everything
> possible, rather than zero.

Maria,
then it's time to debug :-) After the call to
interpolate_boundary_values(), are the values in the 'boundary_values'
map correct? After the call to apply_boundary_values(), are the values
in the 'solution' vector correct? (You can use DataOut to output the
solution vector at this point to look at it graphically.) What are the
values in 'solution' after solving?

If you can't figure out what is wrong, can you shrink your program to
the minimal necessary to illustrate what you are seeing?

Best
W.


--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@colostate.edu
www: http://www.math.colostate.edu/~bangerth/

[Мария Бронзова]

Dear Wolfgang,

Thank you for the valuable hints! I did calculate a very simple system, I applied complex pressures (250,0) on certain boundaries and that is what I get after mapping and after using "apply_boundary_values":

for (auto it = boundary_values.cbegin(); it!= boundary_values.cend(); ++it)

{

std::cout<<it->first<<"  "<<it->second<<std::endl;

}

-------> 

1  (0,0)

2  (0,0)

3 (250,0)

4 (0,0)

5  (0,0)

6  (0,0)

7  (0,0)

8  (250,0)

Which makes perfect sense. But when I print out the solution values with the help of the DataOut for my pressures, I get the following values:

(0,0)

(0,0)

(250,-3.92666e-06)

(0,0)

(0,0)

(0,0)

(0,0)

(250,-5.66168)

I am not quiet sure, what the problem might be. Would you possibly have an idea?

Thank you a lot!

Kind regards,

Mariia

среда, 13 октября 2021 г. в 17:35:12 UTC+2, Wolfgang Bangerth:

[Мария Бронзова]

By the way, the solving procedure is done with the help of the direct solver - SparseDirectUMFPACK.

template <int dim>

  void Poroelasticity<dim>::solve()

  {

         std::cout << "Solving linear system..";

SparseDirectUMFPACK A_direct;

A_direct.initialize(system_matrix);

A_direct.solve(system_matrix, system_rhs);

solution = system_rhs;

time.set_desired_next_step_size(1);

  }

четверг, 14 октября 2021 г. в 11:15:44 UTC+2, Мария Бронзова:

[Wolfgang Bangerth]

On 10/14/21 3:15 AM, Мария Бронзова wrote:
>
> I am not quiet sure, what the problem might be. Would you possibly have
> an idea?

I am not sure either, and can't tell without seeing the code. Like I
said, construct a small testcase where everything extraneous has been
removed (nonlinear loops, time iterations, unnecessary mesh refinement,
etc). Feel free to send such a testcase to the list.

[Мария Бронзова]

Dear Wolfgang,

I believe I answered you privately, but I cannot find my email anywhere. In case you didn"t receive it, I am sending such a test case here again. The testcase is a box porous material model with 8 cells to be able to interpolate BC's and to see the problematics clearlier. The confusing thing might be the Time function, as I am working with the frequency domain, but I simply used it artificially as a counter for frequencies to generate a series of output files. 

The problem appears in the pressure boundary condition, as I am trying to force (250,0) complex-valued boundary values on certain faces. After solving the system, the Dirichlet boundary values are not held, so (250,0) becomes somewhat else.

Please, let me know, if some other information is required!

Thank you a lot!

Kind regards,

Mariia

#include <deal.II/base/quadrature_lib.h>

#include <deal.II/base/logstream.h>

#include <deal.II/base/function.h>

#include <deal.II/lac/block_vector.h>

#include <deal.II/lac/full_matrix.h>

#include <deal.II/lac/block_sparse_matrix.h>

#include <deal.II/lac/solver_cg.h>

#include <deal.II/lac/precondition.h>

#include <deal.II/lac/linear_operator.h>

#include <deal.II/lac/packaged_operation.h>

#include <deal.II/grid/tria.h>

#include <deal.II/grid/grid_generator.h>

#include <deal.II/dofs/dof_handler.h>

#include <deal.II/dofs/dof_renumbering.h>

#include <deal.II/dofs/dof_tools.h>

#include <deal.II/fe/fe_dgq.h>

#include <deal.II/fe/fe_system.h>

#include <deal.II/fe/fe_values.h>

#include <deal.II/numerics/vector_tools.h>

#include <deal.II/numerics/matrix_tools.h>

#include <deal.II/numerics/data_out.h>

#include <deal.II/fe/fe_q.h>

#include <deal.II/lac/sparse_direct.h>

#include <fstream>

#include <iostream>

#include <math.h>

#include <complex>

#include <vector>

#include <deal.II/base/tensor_function.h>

#include <deal.II/base/discrete_time.h>

#include <deal.II/numerics/error_estimator.h>

#include <deal.II/grid/grid_refinement.h>

using namespace std;

namespace Poroacoustics

{

  using namespace dealii;

  template <int dim>

  SymmetricTensor<4, dim, std::complex<double>> get_stress_strain_tensor(const std::complex<double> lambda,

                       const std::complex<double> mu)

  {

  const std::complex<double> nothing(0.0,0.0);

  SymmetricTensor<4, dim, std::complex<double>> stress_strain_tensor;

  for (unsigned int i=0; i< dim; ++i)

  for (unsigned int j=0; j< dim; ++j)

  for (unsigned int k=0; k< dim; ++k)

  for (unsigned int l=0; l< dim; ++l)

  stress_strain_tensor[i][j][k][l]= (((i==k) && (j==l) ? mu : nothing) +

              ((i==l) && (j==k) ? mu : nothing) +

((i==j) && (k==l) ? lambda : nothing));

  return stress_strain_tensor;

  }

  template <int dim>

  inline SymmetricTensor<2, dim> get_strain(const FEValues<dim> &fe_values,

                  const unsigned int shape_func,

const unsigned int q_point)

  {

  SymmetricTensor<2, dim> strain_tensor;

  for (unsigned int i=0; i< dim; ++i)

  strain_tensor[i][i] = fe_values.shape_grad_component(shape_func, q_point, i)[i];

  for (unsigned int i=0; i< dim; ++i)

  for (unsigned int j=i+1; j< dim; ++j)

  strain_tensor[i][j] =

  (fe_values.shape_grad_component(shape_func, q_point, i)[j] +

   fe_values.shape_grad_component(shape_func, q_point, j)[i]) /

   2;

  return strain_tensor;

  }

  template <int dim>

  class Poroelasticity

  {

  public:

    Poroelasticity(const unsigned int degree);

    void run();

  private:

    void make_grid_and_dofs();

    void assemble_system();

    void solve();

    void output_results() const;

    const unsigned int degree;

    Triangulation<dim> triangulation;

    FESystem<dim>      fe;

    DoFHandler<dim>    dof_handler;

    DiscreteTime time;

    BlockSparsityPattern      sparsity_pattern;

    BlockSparseMatrix<std::complex<double>> system_matrix;

    BlockVector<std::complex<double>> solution;

    BlockVector<std::complex<double>> system_rhs;

    static const SymmetricTensor<4, dim> stress_strain_tensor;

static constexpr unsigned int start_freq=100;

static constexpr unsigned int end_freq=101; //[Hz]

  };

  template<int dim>

  class RightHandSide : public Function<dim, std::complex<double>>

  {

  public:

  RightHandSide()

  : Function<dim, std::complex<double>>(dim+1)

  {}

  virtual std::complex<double> value(const Point<dim> &p,

         const unsigned int component = 0) const override;

  };

  template <int dim>

  std::complex<double> RightHandSide<dim>::value(const Point<dim> & p,

               const unsigned int component) const

  {

  const std::complex<double> i = {0.,0.};

  return {0.,0.};

  }

  template <int dim>

  class PressureBoundaryValues : public Function<dim, std::complex<double>>

  {

  public:

  PressureBoundaryValues()

    : Function<dim, std::complex<double>>(dim+1)

      {}

  virtual std::complex<double> value(const Point<dim> &p,

         const unsigned int component = 0) const override;

  };

  template <int dim>

  std::complex<double> PressureBoundaryValues<dim>::value(const Point<dim> & p,

               const unsigned int component) const

  {

  const std::complex<double> i = {0,1};

  return 250.+ 0.0*i;

  }

  template <int dim>

      class DisplBoundaryValues : public Function<dim, std::complex<double>>

      {

      public:

      DisplBoundaryValues()

        : Function<dim, std::complex<double>>(dim+1)

          {}

      virtual std::complex<double> value(const Point<dim> &p,

               const unsigned int component = 0) const override;

      virtual void vector_value(const Point<dim> &p,

      Vector<std::complex<double>> &  value) const override;

      };

    template <int dim>

  std::complex<double> DisplBoundaryValues<dim>::value(const Point<dim> & p,

     const unsigned int component) const

{

   //const std::complex<double> i = {0,1};

   return {0.0,0.0};

}

    template <int dim>

    void DisplBoundaryValues<dim>::vector_value(const Point<dim> &p,

                                           Vector<std::complex<double>> & values) const

    {

      const std::complex<double> i = {0,1};

      values(0) = DisplBoundaryValues<dim>::value(p, 0);

      values(1) = DisplBoundaryValues<dim>::value(p, 1);

      values(2) = 0.0+i*0.0;

    }

  template <int dim>

  Poroelasticity<dim>::Poroelasticity(const unsigned int degree)

    : degree(degree)

  , triangulation(Triangulation<dim>::maximum_smoothing)

    , fe(FE_Q<dim>(degree+1),dim, FE_Q<dim>(degree),1)

    , dof_handler(triangulation)

    , time(/*start time*/ 0., /*end time*/ 1.) //Here the second argument needs to be manually adjusted

  {}

  template <int dim>

  void Poroelasticity<dim>::make_grid_and_dofs()

  {

std::vector<unsigned int> subdivisions(dim,1);

subdivisions[0]=2;

subdivisions[1]=2;

subdivisions[2]=2;

const Point<dim> bottom_left = (dim == 2?

Point<dim>(-1,-1) :

Point<dim>(-0.5, 0, -0.5));

const Point<dim> top_right = (dim == 2?

Point<dim>(1, 0) :

Point<dim>(0.5, 1, 0));

GridGenerator::subdivided_hyper_rectangle(triangulation,

   subdivisions,

   bottom_left,

   top_right,

   false);

for (const auto &cell : triangulation.active_cell_iterators())

for (const auto &face : cell->face_iterators())

if ((face->center()[dim-1] == 0))

face->set_all_boundary_ids(1);

else if  (face->center()[0]==-0.5)

face->set_all_boundary_ids(2);

else if (face->boundary_id()!=2 && face->boundary_id()!=1

&&face->at_boundary())

{

face->set_all_boundary_ids(3);

}

    dof_handler.distribute_dofs(fe);

    DoFRenumbering::Cuthill_McKee(dof_handler);

    std::vector<unsigned int> block_component(dim+1, 0);

    block_component[dim]=1;

    DoFRenumbering::component_wise(dof_handler, block_component);

    const std::vector<types::global_dof_index> dofs_per_block =

          DoFTools::count_dofs_per_fe_block(dof_handler, block_component);

    const unsigned int n_u = dofs_per_block[0],

                       n_p = dofs_per_block[1];

    std::cout << "Number of active cells: " << triangulation.n_active_cells()

              << std::endl

              << "Total number of cells: " << triangulation.n_cells()

              << std::endl

              << "Number of degrees of freedom: " << dof_handler.n_dofs()

              << " (" << n_u << '+' << n_p << ')' << std::endl;

    BlockDynamicSparsityPattern dsp(2, 2);

    dsp.block(0, 0).reinit(n_u, n_u);

    dsp.block(1, 0).reinit(n_p, n_u);

    dsp.block(0, 1).reinit(n_u, n_p);

    dsp.block(1, 1).reinit(n_p, n_p);

    dsp.collect_sizes();

    DoFTools::make_sparsity_pattern(dof_handler, dsp); //, constraints, false);

    sparsity_pattern.copy_from(dsp);

    system_matrix.reinit(sparsity_pattern);

    solution.reinit(2);

    solution.block(0).reinit(n_u);

    solution.block(1).reinit(n_p);

    solution.collect_sizes();

    system_rhs.reinit(2);

    system_rhs.block(0).reinit(n_u);

    system_rhs.block(1).reinit(n_p);

    system_rhs.collect_sizes();

  }

  template <int dim>

  void Poroelasticity<dim>::assemble_system()

  {

    QGauss<dim>     quadrature_formula(degree + 2);

    QGauss<dim - 1> face_quadrature_formula(degree + 2);

    FEValues<dim>     fe_values(fe,

                            quadrature_formula,

                            update_values | update_gradients |

                              update_quadrature_points | update_JxW_values);

    FEFaceValues<dim> fe_face_values(fe,

                                     face_quadrature_formula,

                                     update_values | update_normal_vectors |

                                       update_quadrature_points |

                                       update_JxW_values);

    const unsigned int dofs_per_cell   = fe.n_dofs_per_cell();

    const unsigned int n_q_points      = quadrature_formula.size();

    std::cout <<"N_q_points="<< n_q_points<< std::endl;

    const unsigned int n_face_q_points = face_quadrature_formula.size();

    std::cout <<"N_face_q_points="<< n_face_q_points<< std::endl;

    constexpr double porosity = 0.96;

constexpr double resistivity  = 87000; //[N*s/m^4]

constexpr double tortuosity  = 2.52;

constexpr double visc_length  = 3.7e-05; //[m]

constexpr double therm_length  = 0.000116; //[m]

constexpr double loss_factor  = 0.01;

constexpr double E_modulus = 56000000; //[Pa]

constexpr double poisson = 0.3;

constexpr double frame_density = 31; //[kg/m^3]

constexpr double fluid_density = 1.2; //[kg/m^3]

constexpr double Prandtl_number = 0.71; //[-]

const std::complex<double> i = {0,1};

const std::complex<double> shear_modulus = (220.+i*22.)*1e-4; //[N*m-2]

//constexpr double thermal_conductivity = 2.6*1e-2; //[w*m-1*k-1]

constexpr double specific_heat_ratio = 1.4; //[cp/cv]

constexpr double ambient_pressure = 1.0132*1e5; //[Pa]

std::vector<int> freq(end_freq-0);

std::vector<double> omega(freq.size());

std::vector<std::complex<double>> visc_factor(freq.size());

std::vector<std::complex<double>> dyn_tortuosity(freq.size());

std::vector<std::complex<double>> density11(freq.size());

std::vector<std::complex<double>> density12(freq.size());

std::vector<std::complex<double>> density22(freq.size());

std::vector<std::complex<double>> density0(freq.size());

std::vector<std::complex<double>> bulk_fluid(freq.size());

std::vector<std::complex<double>> A(freq.size());

std::vector<std::complex<double>> Q(freq.size());

std::vector<std::complex<double>> R(freq.size());

std::vector<std::complex<double>> gamma(freq.size());

const std::complex<double> bulk_sk_vacuum = 2.*shear_modulus*(poisson+1.)/(3.*(1.-2.*poisson)); //[]

const std::complex<double> bulk_elast_solid = poisson*E_modulus/((poisson+1.)*(1.-2.*poisson))+2.*shear_modulus/3.; //[]

FullMatrix<std::complex<double>> local_matrix(dofs_per_cell, dofs_per_cell);

Vector<std::complex<double>>     local_rhs(dofs_per_cell);

for (unsigned int j=0; j<end_freq-start_freq; ++j)

{

freq[j]=j+start_freq;

//std::cout<< freq[j]<<std::endl;

omega[j]=2*M_PI*freq[j];

//std::cout<< omega[j]<<std::endl;

visc_factor[j]=std::sqrt(1.+i*4.*std::pow(tortuosity, 2)*loss_factor*fluid_density*omega[j]/

(std::pow(resistivity*visc_length*porosity, 2)));

dyn_tortuosity[j]=tortuosity-i*porosity*resistivity*visc_factor[j]/(fluid_density*omega[j]);

//std::cout<< dyn_tortuosity[j]<<std::endl;

density11[j]=(1.-porosity)*frame_density+porosity*fluid_density*(dyn_tortuosity[j]-1.);

density12[j]=-porosity*fluid_density*(dyn_tortuosity[j]-1.);

density22[j]=porosity*fluid_density*dyn_tortuosity[j];

density0[j]=density11[j]-std::pow(density12[j],2)/density22[j];

//std::cout<< density0[j]<<std::endl;

bulk_fluid[j]=ambient_pressure/(1.-(specific_heat_ratio-1.)/(specific_heat_ratio*(8.*loss_factor*

std::sqrt(1.+std::pow(therm_length/4.,2)*i*omega[j]*fluid_density*std::pow(Prandtl_number,2)/loss_factor)/

(fluid_density*std::pow(Prandtl_number*therm_length,2)*i*omega[j])+1.)));

//std::cout<< bulk_fluid[j]<<std::endl;

A[j]=((1-porosity)*(1-porosity-bulk_sk_vacuum/bulk_elast_solid)*

bulk_elast_solid+porosity*(bulk_elast_solid/bulk_fluid[j])*bulk_sk_vacuum)/

(1.-porosity-bulk_sk_vacuum/bulk_elast_solid+porosity*bulk_elast_solid/bulk_fluid[j])-2.*shear_modulus/3.;

//std::cout<< A[j]<<std::endl;

Q[j]=((1-porosity-bulk_sk_vacuum/bulk_elast_solid)*porosity*bulk_elast_solid)/

(1-porosity-bulk_sk_vacuum/bulk_elast_solid+porosity*bulk_elast_solid/bulk_fluid[j]);

//std::cout<< Q[j]<<std::endl;

R[j]=std::pow(porosity,2)*bulk_elast_solid/

(1-porosity-bulk_sk_vacuum/bulk_elast_solid+porosity*bulk_elast_solid/bulk_fluid[j]);

//std::cout<< R[j]<<std::endl;

gamma[j]=porosity*(density12[j]/density22[j] - Q[j]/R[j]);

//std::cout<< gamma[j]<<std::endl;

}

    std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);

    const RightHandSide<dim> right_hand_side;

    const PressureBoundaryValues<dim> pressure_boundary_values;

    const DisplBoundaryValues<dim> displ_boundary_values;

    std::vector<std::complex<double>>         rhs_values(n_q_points);

    std::vector<std::complex<double>>         p_boundary_values(n_face_q_points);

    std::vector<Vector<std::complex<double>>> d_boundary_values(n_face_q_points, Vector<std::complex<double>>(dim + 1));

    const FEValuesExtractors::Vector velocities(0);

    const FEValuesExtractors::Scalar pressure(dim);

 for (unsigned int k=0; k< end_freq-start_freq; ++k)

    {

for (const auto &cell : dof_handler.active_cell_iterators())

  {

fe_values.reinit(cell);

local_matrix = 0;

local_rhs    = 0;

right_hand_side.value_list(fe_values.get_quadrature_points(),rhs_values);

for (unsigned int q = 0; q < n_q_points; ++q)

  for (unsigned int i = 0; i < dofs_per_cell; ++i)

{

  const Tensor<1, dim> phi_i_u = fe_values[velocities].value(i, q);

// const double div_phi_i_u = fe_values[velocities].divergence(i, q);

  const double phi_i_p     = fe_values[pressure].value(i, q);

  const Tensor<1,dim> div_phi_i_p = fe_values[pressure].gradient(i, q);

  for (unsigned int j = 0; j < dofs_per_cell; ++j)

{

  const Tensor<1, dim> phi_j_u =

fe_values[velocities].value(j, q);

// const double div_phi_j_u =

// fe_values[velocities].divergence(j, q);

  const double phi_j_p = fe_values[pressure].value(j, q);

  const Tensor<1,dim> div_phi_j_p = fe_values[pressure].gradient(j, q);

  const SymmetricTensor<2, dim>

    eps_phi_i = get_strain(fe_values, i, q),

    eps_phi_j = get_strain(fe_values, j, q);

  const SymmetricTensor<4, dim, std::complex<double>>

    material_matrix = get_stress_strain_tensor<dim>(/*lambda = */ A[k]-std::pow(Q[k],2)/R[k],

                        /*mu = */ 2.*shear_modulus);

  local_matrix(i, j) +=

(eps_phi_i * material_matrix * eps_phi_j // transpose

- std::pow(omega[k],2)*density0[k]*phi_i_u*phi_j_u  //

- gamma[k]*div_phi_i_p*phi_j_u

+ std::pow(porosity,2)*div_phi_i_p*div_phi_j_p/(density22[k]*std::pow(omega[k],2))

- std::pow(porosity,2)*phi_i_p*phi_j_p/R[k]

- gamma[k]*div_phi_i_p*phi_j_u)                //

* fe_values.JxW(q);

}

  local_rhs(i) += rhs_values[q]* fe_values.JxW(q);

}

for (const auto &face : cell->face_iterators())

  if (face->boundary_id()==1)

  {

  //std::cout<<"Boundary ID = 1 detected"<<std::endl;

  fe_face_values.reinit(cell, face);

  pressure_boundary_values.value_list(

  fe_face_values.get_quadrature_points(), p_boundary_values);

  displ_boundary_values.vector_value_list(fe_face_values.get_quadrature_points(),

      d_boundary_values);

  for (unsigned int q = 0; q < n_face_q_points; ++q)

  for (unsigned int i = 0; i < dofs_per_cell; ++i)

  {

  local_rhs(i) += (1.-porosity*(1.+Q[k]/R[k]))

  *p_boundary_values[q]

  *fe_face_values[velocities].value(i,q)

  *fe_face_values.normal_vector(q)

  *fe_face_values.JxW(q);

  }

  }

  else if (face->boundary_id()==2)

  {

  std::cout<<"Boundary ID = 2 detected"<<std::endl;

  fe_face_values.reinit(cell, face);

  displ_boundary_values.vector_value_list(fe_face_values.get_quadrature_points(),

    d_boundary_values);

  for (unsigned int q = 0; q < n_face_q_points; ++q)

  for (unsigned int i = 0; i < dofs_per_cell; ++i)

  for (unsigned int c = 0; c < dim; ++c)

  {

          const unsigned int component_i = fe.system_to_component_index(i).first;

  local_rhs(i) += -porosity*(1.+Q[k]/R[k])

  *d_boundary_values[q](component_i)

  *fe_face_values[pressure].value(i,q)

  *fe_face_values.JxW(q);

  }

  }

cell->get_dof_indices(local_dof_indices);

for (unsigned int i = 0; i < dofs_per_cell; ++i)

  for (unsigned int j = 0; j < dofs_per_cell; ++j)

system_matrix.add(local_dof_indices[i],

  local_dof_indices[j],

  local_matrix(i, j));

for (unsigned int i = 0; i < dofs_per_cell; ++i)

  system_rhs(local_dof_indices[i]) += local_rhs(i);

    }//active cell iterator loop

FEValuesExtractors::Scalar pressure(dim);

FEValuesExtractors::Vector velocities(0);

std::map<types::global_dof_index, std::complex<double>> boundary_values;

VectorTools::interpolate_boundary_values(dof_handler,

3,

RightHandSide<dim>(),

boundary_values,

fe.component_mask(pressure));

VectorTools::interpolate_boundary_values(dof_handler,

1,

PressureBoundaryValues<dim>(),

boundary_values,

fe.component_mask(pressure));

VectorTools::interpolate_boundary_values(dof_handler,

2,

DisplBoundaryValues<dim>(),

boundary_values,

fe.component_mask(velocities));

for (auto it = boundary_values.cbegin(); it!= boundary_values.cend(); ++it)

{

std::cout<<it->first<<"  "<<it->second<<std::endl;

}

    MatrixTools::apply_boundary_values(boundary_values,

   system_matrix,

   solution,

   system_rhs);

//     std::cout<<"--------------------------------"<<std::endl;

//     for (auto it = boundary_values.cbegin(); it!= boundary_values.cend(); ++it)

// {

// std::cout<<it->first<<"  "<<it->second<<std::endl;

// }

std::cout << "Frequencystep " << time.get_step_number() + 1 << std::endl;

std::ofstream out("out-" +

                             Utilities::int_to_string(time.get_step_number(), 4) +

                             ".txt");

system_matrix.print_formatted(out,3,true,0,"(0,0)");

//std::ofstream out1("out1.txt");

std::ofstream out1("out1-" +

                             Utilities::int_to_string(time.get_step_number(), 4) +

                             ".txt");

system_rhs.print(out1,3,true,true);

if (time.is_at_end() != true)

{

solve();

output_results();

time.advance_time();

}

  }//frequency loop

  }

  template <int dim>

  void Poroelasticity<dim>::solve()

  {

  std::cout << "Solving linear system..";

  /*for (unsigned int k=0; k<end_freq-start_freq; ++k)

  {*/

SparseDirectUMFPACK A_direct;

A_direct.initialize(system_matrix);

//A_direct.vmult(solution, system_rhs);

A_direct.solve(system_matrix, system_rhs);

solution = system_rhs;

time.set_desired_next_step_size(1);

  //}

  }

  namespace DataPostprocessors

    {

      template <int dim>

      class ComplexAmplitude : public DataPostprocessorScalar<dim>

      {

      public:

        ComplexAmplitude();

        virtual void evaluate_vector_field(

          const DataPostprocessorInputs::Vector<dim> &inputs,

          std::vector<Vector<double>> &computed_quantities) const override;

      };

      template <int dim>

      ComplexAmplitude<dim>::ComplexAmplitude()

        : DataPostprocessorScalar<dim>("Amplitude", update_values)

      {}

      template <int dim>

      void ComplexAmplitude<dim>::evaluate_vector_field(

        const DataPostprocessorInputs::Vector<dim> &inputs,

        std::vector<Vector<double>> &computed_quantities) const

      {

        Assert(computed_quantities.size() == inputs.solution_values.size(),

               ExcDimensionMismatch(computed_quantities.size(),

                                    inputs.solution_values.size()));

        for (unsigned int q = 0; q < computed_quantities.size(); ++q)

          {

            Assert(computed_quantities[q].size() == 1,

                   ExcDimensionMismatch(computed_quantities[q].size(), 1));

            Assert(inputs.solution_values[q].size() == 8,

                   ExcDimensionMismatch(inputs.solution_values[q].size(), 8));

            const std::complex<double> psi(inputs.solution_values[q](3),

                                           inputs.solution_values[q](7));

            std::cout<<psi<<std::endl;

            computed_quantities[q](0) = std::norm(psi);

          }

      }

      } // namespace DataPostprocessors

  template <int dim>

  void Poroelasticity<dim>::output_results() const

  {

    std::vector<std::string> solution_names(1, "u1");

    //std::vector<std::string> solution_names(dim, "u");

    solution_names.emplace_back("u2");

    solution_names.emplace_back("u3");

    solution_names.emplace_back("p");

    std::vector<DataComponentInterpretation::DataComponentInterpretation>

      interpretation(dim,

                     //DataComponentInterpretation::component_is_part_of_vector);

        DataComponentInterpretation::component_is_scalar);

    interpretation.push_back(DataComponentInterpretation::component_is_scalar);

    const DataPostprocessors::ComplexAmplitude<dim> complex_magnitude;

    DataOut<dim> data_out;

    data_out.attach_dof_handler(dof_handler);

    data_out.add_data_vector(dof_handler,

                             solution,

                             solution_names,

                             interpretation);

    data_out.add_data_vector(solution, complex_magnitude);

    data_out.build_patches();

    //std::ofstream output("solution.vtu");

    std::ofstream output("solution-" +

                             Utilities::int_to_string(time.get_step_number(), 4) +

                             ".vtu");

    data_out.write_vtu(output);

  }

  template <int dim>

  void Poroelasticity<dim>::run()

  {

    make_grid_and_dofs();

assemble_system();

  }

} // namespace Poroacoustics

int main()

{

  #ifdef DEBUG

{

std::cout<<"HEllo";

}

  #endif

  try

    {

      using namespace Poroacoustics;

      const unsigned int     fe_degree = 1;

      Poroelasticity<3> poroelasticity(fe_degree);

      poroelasticity.run();

    }

  catch (std::exception &exc)

    {

      std::cerr << std::endl

                << std::endl

                << "----------------------------------------------------"

                << std::endl;

      std::cerr << "Exception on processing: " << std::endl

                << exc.what() << std::endl

                << "Aborting!" << std::endl

                << "----------------------------------------------------"

                << std::endl;

      return 1;

    }

  catch (...)

    {

      std::cerr << std::endl

                << std::endl

                << "----------------------------------------------------"

                << std::endl;

      std::cerr << "Unknown exception!" << std::endl

                << "Aborting!" << std::endl

                << "----------------------------------------------------"

                << std::endl;

      return 1;

    }

  return 0;

}

пятница, 15 октября 2021 г. в 01:06:59 UTC+2, Wolfgang Bangerth:

[Wolfgang Bangerth]

On 10/21/21 7:33 AM, Мария Бронзова wrote:
>
> I believe I answered you privately, but I cannot find my email anywhere.
> In case you didn"t receive it, I am sending such a test case here again.
> The testcase is a box porous material model with 8 cells to be able to
> interpolate BC's and to see the problematics clearlier. The confusing
> thing might be the Time function, as I am working with the frequency
> domain, but I simply used it artificially as a counter for frequencies
> to generate a series of output files.
>
> The problem appears in the pressure boundary condition, as I am trying
> to force (250,0) complex-valued boundary values on certain faces. After
> solving the system, the Dirichlet boundary values are not held, so
> (250,0) becomes somewhat else.
>
> Please, let me know, if some other information is required!

Maria -- the program is still too large. For debugging these kinds of
problems, it is often useful to strip everything that is not necessary
out. Set all of your coefficients to one (i.e., omit them from the code
altogether), remove the loop over frequencies, remove the postprocessing
steps other than graphical output, etc. I bet that you can reduce your
current 724 lines to 250 this way, and it will be much easier to see
what could possibly go wrong in your code then.

[Мария Бронзова]

Dear Wolfgang,

I am sorry for not having answered to your reply. I started using the constraints at some point and they worked well for my problem. But I came back to analyzing the issue again now.

I eliminated all the coefficients, unnecessary loops and so on, as you suggested. The code has still somewhat 450 lines, but I tried to leave only the core parts. As I was eliminating things, I found out, where the problem occurs, but I have not an explanation for it. 

Basically, now I only have one additional factor in my local_matrix (apart from it's "skeleton") within the Assembly module - density22, which is a complex value (0,1). When I remove it, the Dirichlet BC get applied as expected through the MatrixTools::apply_boundary_values function. But when i leave it in the system - boundary conditions are no longer correct after solving the system.

Would you please give a hint, what is going on here?

Thank you for your support!

Here is the simplified code. 

Kind regards,

Mariia

    

понедельник, 25 октября 2021 г. в 22:21:31 UTC+2, Wolfgang Bangerth:

[Wolfgang Bangerth]

On 12/2/21 7:10 AM, Мария Бронзова wrote:
>
> Basically, now I only have one additional factor in my local_matrix (apart
> from it's "skeleton") within the Assembly module - density22, which is a
> complex value (0,1). When I remove it, the Dirichlet BC get applied as
> expected through the MatrixTools::apply_boundary_values function. But when i
> leave it in the system - boundary conditions are no longer correct after
> solving the system.
>
> Would you please give a hint, what is going on here?

I *think* what might be happening is that you are running into a bug that
someone else just a few weeks ago discovered:
https://github.com/dealii/dealii/issues/12916
The patch for it is here:
https://github.com/dealii/dealii/pull/13041

I realize that that must have been very frustrating to debug :-( Can you try
and see whether the issue you see goes away if you apply that patch and
recompile deal.II with it?

Best
Wolfgang

[Мария Бронзова]

Dear Wolfgang, 

it was indeed frustrating, took me a while to detect the section, where it starts to break down. But it is what it is.

But in the meantime i recompiled the library and it looks good now, what a relief! Compared the results to the reference I have and it's matching now.

Thank you a lot for the solution, good that this forum exists!

Kind regards,

Mariia

четверг, 9 декабря 2021 г. в 04:44:01 UTC+1, Wolfgang Bangerth: