Jaekwang Kim
template <int dim>
class RightHandSide : public Function<dim>
{
public:
RightHandSide () : Function<dim>(dim+1) {}
virtual void vector_value (const Point<dim> &p,
Vector<double> &value) const;
};
template <int dim>
void
RightHandSide<dim>::vector_value (const Point<dim> &p,
Vector<double> &values) const
{
double n = 0.5 ;
double x =p[0]; double y=p[1];
double f_x=2*x + 7*y - 4*Power(E,3*(Power(-0.5 + x,2) + Power(-0.5 + y,2)))*(-0.5 + n)*
(1 + 2*Power(-0.5 + x,2))*(-0.5 + x)*
(240. - 320.*x + 448.*Power(x,2) - 256.*Power(x,3) + 128.*Power(x,4) - 192.*y +
320.*Power(y,2) - 256.*Power(y,3) + 128.*Power(y,4))*
Power(1 + (Power(E,2*(Power(-0.5 + x,2) + Power(-0.5 + y,2)))*
(64.*Power(0.75 - 1.*x + Power(x,2),2) + 64.*Power(0.75 - 1.*y + Power(y,2),2)))
/2.,-1.5 + n) - 16.*Power(E,Power(-0.5 + x,2) + Power(-0.5 + y,2))*(-0.5 + x)*
Power(1 + (Power(E,2*(Power(-0.5 + x,2) + Power(-0.5 + y,2)))*
(64.*Power(0.75 - 1.*x + Power(x,2),2) + 64.*Power(0.75 - 1.*y + Power(y,2),2)))
/2.,-0.5 + n) - 8*Power(E,Power(-0.5 + x,2) + Power(-0.5 + y,2))*
(1 + 2*Power(-0.5 + x,2))*(-0.5 + x)*
Power(1 + (Power(E,2*(Power(-0.5 + x,2) + Power(-0.5 + y,2)))*
(64.*Power(0.75 - 1.*x + Power(x,2),2) + 64.*Power(0.75 - 1.*y + Power(y,2),2)))
/2.,-0.5 + n);
double f_y=7*x + 8.*Power(E,3*(Power(-0.5 + x,2) + Power(-0.5 + y,2)))*(-0.5 + n)*(-0.5 + y)*
(0.75 - 1.*y + Power(y,2))*(240. - 192.*x + 320.*Power(x,2) - 256.*Power(x,3) +
128.*Power(x,4) - 320.*y + 448.*Power(y,2) - 256.*Power(y,3) + 128.*Power(y,4))*
Power(1 + (Power(E,2*(Power(-0.5 + x,2) + Power(-0.5 + y,2)))*
(64.*Power(0.75 - 1.*x + Power(x,2),2) + 64.*Power(0.75 - 1.*y + Power(y,2),2)))
/2.,-1.5 + n) + 16.*Power(E,Power(-0.5 + x,2) + Power(-0.5 + y,2))*(-0.5 + y)*
Power(1 + (Power(E,2*(Power(-0.5 + x,2) + Power(-0.5 + y,2)))*
(64.*Power(0.75 - 1.*x + Power(x,2),2) + 64.*Power(0.75 - 1.*y + Power(y,2),2)))
/2.,-0.5 + n) + 16.*Power(E,Power(-0.5 + x,2) + Power(-0.5 + y,2))*(-0.5 + y)*
(0.75 - 1.*y + Power(y,2))*Power(1 +
(Power(E,2*(Power(-0.5 + x,2) + Power(-0.5 + y,2)))*
(64.*Power(0.75 - 1.*x + Power(x,2),2) + 64.*Power(0.75 - 1.*y + Power(y,2),2)))
/2.,-0.5 + n);
values(0) = f_x;
values(1) = f_y;
values(2) = 0;
}
template <int dim> // Described in Cartesian Coordinate
void ShearProblem<dim>::assemble_system ()
{
system_matrix=0;
system_rhs=0;
// Variable Coefficient
double shear_rate;
double viscosity;
const MappingQ<dim> mapping (degree);
QGauss<dim> quadrature_formula(4*degree+1);
FEValues<dim> fe_values (mapping, fe, quadrature_formula,
update_values |
update_quadrature_points |
update_JxW_values |
update_gradients);
const unsigned int dofs_per_cell = fe.dofs_per_cell;
const unsigned int n_q_points = quadrature_formula.size();
FullMatrix<double> local_matrix (dofs_per_cell, dofs_per_cell);
Vector<double> local_rhs (dofs_per_cell);
std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
const RightHandSide<dim> right_hand_side;
std::vector<Vector<double> > rhs_values (n_q_points, Vector<double>(dim+1));
const FEValuesExtractors::Vector velocities (0);
const FEValuesExtractors::Scalar xvel (0);
const FEValuesExtractors::Scalar pressure (dim);
std::vector<SymmetricTensor<2,dim> > symgrad_phi_u (dofs_per_cell);
std::vector<double> div_phi_u (dofs_per_cell);
std::vector<double> phi_p (dofs_per_cell);
std::vector<SymmetricTensor<2,dim> > local_previous_symgrad_phi_u (n_q_points);
std::vector<double> local_previous_solution (n_q_points);
typename DoFHandler<dim>::active_cell_iterator
cell = dof_handler.begin_active(),
endc = dof_handler.end();
for (; cell!=endc; ++cell)
{
fe_values.reinit (cell);
local_matrix = 0;
local_rhs = 0;
right_hand_side.vector_value_list(fe_values.get_quadrature_points(),rhs_values);
fe_values[velocities].get_function_symmetric_gradients (previous_solution, local_previous_symgrad_phi_u);
for (unsigned int q=0; q<n_q_points; ++q)
{
// AT Each Quadrature Point, Viscosity is calcualted from previous solution
shear_rate = std::sqrt( local_previous_symgrad_phi_u[q]* local_previous_symgrad_phi_u[q] *2 );
//HB model
viscosity = std::pow( 1+pow(shear_rate,2) , (power_n-1)/2);
for (unsigned int k=0; k<dofs_per_cell; ++k)
{
symgrad_phi_u[k] = fe_values[velocities].symmetric_gradient (k, q);
div_phi_u[k] = fe_values[velocities].divergence (k, q);
phi_p[k] = fe_values[pressure].value (k, q);
}
for (unsigned int i=0; i<dofs_per_cell; ++i)
{
for (unsigned int j=0; j<dofs_per_cell; ++j)
{
local_matrix(i,j) += (2 * viscosity * (symgrad_phi_u[i] * symgrad_phi_u[j])
- div_phi_u[i] * phi_p[j]
- phi_p[i] * div_phi_u[j]
+ phi_p[i] * phi_p[j]) * fe_values.JxW(q);
}
}
for (unsigned int i=0; i<dofs_per_cell; ++i)
{
local_rhs(i) += fe_values.shape_value(i,q) *
rhs_values[q] *
fe_values.JxW(q);
}
}
cell->get_dof_indices (local_dof_indices);
constraints.distribute_local_to_global (local_matrix, local_rhs,
local_dof_indices,
system_matrix, system_rhs);
}
/Users/jaekwangjk/repos/simpleshear/shear.cc:609:60: error: invalid operands to
binary expression ('double' and 'value_type' (aka
'dealii::Vector<double>'))
local_rhs(i) += fe_values.shape_value(i,q) *
~~~~~~~~~~~~~~~~~~~~~~~~~~ ^
/Users/jaekwangjk/repos/simpleshear/shear.cc:879:9: note: in instantiation of
member function 'MyStokes::ShearProblem<2>::assemble_system' requested
here
assemble_system ();
^
/Users/jaekwangjk/repos/simpleshear/shear.cc:1054:22: note: in instantiation of
member function 'MyStokes::ShearProblem<2>::run' requested here
simple_shear.run ();
^
/Users/jaekwangjk/Programs/dealii/include/deal.II/base/complex_overloads.h:39:1: note:
candidate template ignored: could not match 'complex<type-parameter-0-0>'
against 'const double'
operator*(const std::complex<T> &left, const std::complex<U> &right)
^
/Users/jaekwangjk/Programs/dealii/include/deal.II/base/complex_overloads.h:57:1: note:
candidate template ignored: could not match 'complex<type-parameter-0-0>'
against 'const double'
operator*(const std::complex<T> &left, const U &right)
^
/Users/jaekwangjk/Programs/dealii/include/deal.II/base/complex_overloads.h:75:1: note:
candidate template ignored: could not match 'complex' against 'Vector'
operator*(const T &left, const std::complex<U> &right)
^
/Users/jaekwangjk/Programs/dealii/include/deal.II/base/tensor.h:1232:1: note:
candidate template ignored: could not match 'Tensor' against 'Vector'
operator * (const Other object,
^
/Users/jaekwangjk/Programs/dealii/include/deal.II/base/tensor.h:1250:1: note:
candidate template ignored: could not match 'Tensor<0, dim,
type-parameter-0-1>' against 'const double'
operator * (const Tensor<0,dim,Number> &t,
^
/Users/jaekwangjk/Programs/dealii/include/deal.II/base/tensor.h:1269:1: note:
candidate template ignored: could not match 'Tensor<0, dim,
type-parameter-0-1>' against 'const double'
operator * (const Tensor<0, dim, Number> &src1,
^
/Users/jaekwangjk/Programs/dealii/include/deal.II/base/tensor.h:1335:1: note:
candidate template ignored: could not match 'Tensor<rank, dim,
type-parameter-0-2>' against 'const double'
operator * (const Tensor<rank,dim,Number> &t,
^
/Users/jaekwangjk/Programs/dealii/include/deal.II/base/tensor.h:1361:1: note:
candidate template ignored: could not match 'Tensor' against 'Vector'
operator * (const Number factor,
^
/Users/jaekwangjk/Programs/dealii/include/deal.II/base/tensor.h:1468:1: note:
candidate template ignored: could not match 'Tensor<rank_1, dim,
type-parameter-0-3>' against 'const double'
operator * (const Tensor<rank_1, dim, Number> &src1,
^
/Users/jaekwangjk/Programs/dealii/include/deal.II/base/point.h:482:1: note:
candidate template ignored: could not match 'Point' against 'Vector'
operator * (const OtherNumber factor,
^
/Users/jaekwangjk/Programs/dealii/include/deal.II/base/symmetric_tensor.h:2635:1: note:
candidate template ignored: could not match 'SymmetricTensor<rank, dim,
type-parameter-0-2>' against 'const double'
operator * (const SymmetricTensor<rank,dim,Number> &t,
^
/Users/jaekwangjk/Programs/dealii/include/deal.II/base/symmetric_tensor.h:2655:1: note:
candidate template ignored: could not match 'SymmetricTensor' against
'Vector'
operator * (const Number factor,
^
/Users/jaekwangjk/Programs/dealii/include/deal.II/base/symmetric_tensor.h:2709:1: note:
candidate template ignored: could not match 'SymmetricTensor<rank, dim,
type-parameter-0-2>' against 'const double'
operator * (const SymmetricTensor<rank,dim,Number> &t,
^
/Users/jaekwangjk/Programs/dealii/include/deal.II/base/symmetric_tensor.h:2736:1: note:
candidate template ignored: could not match 'SymmetricTensor' against
'Vector'
operator * (const Number factor,
^
/Users/jaekwangjk/Programs/dealii/include/deal.II/base/symmetric_tensor.h:2772:1: note:
candidate template ignored: could not match 'SymmetricTensor<rank, dim,
double>' against 'const double'
operator * (const SymmetricTensor<rank,dim> &t,
^
/Users/jaekwangjk/Programs/dealii/include/deal.II/base/symmetric_tensor.h:2791:1: note:
candidate template ignored: could not match 'SymmetricTensor' against
'Vector'
operator * (const double factor,
^
/Users/jaekwangjk/Programs/dealii/include/deal.II/base/symmetric_tensor.h:3086:1: note:
candidate template ignored: could not match 'SymmetricTensor<2, dim,
type-parameter-0-1>' against 'const double'
operator * (const SymmetricTensor<2,dim,Number> &src1,
^
1 error generated.
Wolfgang Bangerth
> */_local_rhs(i) += fe_values.shape_value(i,q) *_/*
>
> */_ rhs_values[q] *_/*
>
> */_ fe_values.JxW(q);_/*
>
rhs_values[q] is a Vector<double>. You need to say which component of it
you want to multiply with. (Hint: The correct component is
fe.system_to_component_index(i).first).
Best
W.
--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@colostate.edu
www: http://www.math.colostate.edu/~bangerth/
Jaekwang Kim
2017년 2월 8일 수요일 오후 1시 40분 46초 UTC-6, Wolfgang Bangerth 님의 말:
Wolfgang Bangerth
>
> I am using deal-ii to analyze Non-newtonian fluids, which varies viscosity.
> For example, my viscosity eta is function of (u) field. I am solving
> this with iteration method.. (Picard)
> After coding, I am testing my code with method of manufacture.
>
> For squared domain, I implemented all boundary condition with my assumed
> solution and impose source term also to numerically generate my assumed
> solution
> however, the problem I got is.... also follow
>
>
> While my velocity field is qualitatively similar to manufactured solution,
> the pressure field is not following the exact form as seen below figure.
> It is very strange that......I can approach u-soltution form while my
> pressure does not....
>
>
> In addition , any mesh refinement would not make more accurate solution
> now.
>
> If you were me, what would you check more to find out what is problem?
have a constant viscosity? If you iterate, do you get the solution after
one iteration?
I will also note that your manufactured solution is symmetric, but your
computational solution is not. This provides you with a powerful way
because you don't actually need to look at the convergence, it's enough
to check that on a *coarse* mesh the solution is *symmetric*. It may be
that solution is already non-symmetric on 1 cell, or on a 2x2 mesh -- in
which case you already know that something is wrong. The question then
is why that is so -- is your rhs correct, for example?
Jaekwang Kim
I would make the problem simpler. Do you get the right solution if you
have a constant viscosity? If you iterate, do you get the solution after
one iteration?
I will also note that your manufactured solution is symmetric, but your
computational solution is not. This provides you with a powerful way
because you don't actually need to look at the convergence, it's enough
to check that on a *coarse* mesh the solution is *symmetric*. It may be
that solution is already non-symmetric on 1 cell, or on a 2x2 mesh -- in
which case you already know that something is wrong. The question then
is why that is so -- is your rhs correct, for example?
<2x2 mesh>
<4x4 mesh>
Jaekwang Kim
Wolfgang Bangerth
>
> The manufactured solution of velocities has not satisfies the continuity
> equation.
> (i.e. manufactured solution does not satisfy, div.u=0)...