Hi all,
I was testing my code, with method of manufactured solution.
I have a source term, which is vector value, defined as ..
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;
}
and I wanted to assemble RHS in assemble_system function
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);
}
but I have compile error when I tried to run the code.
which is like ...
*/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.
I would appreciate if someone tells me what I am doing wrong...
Thank you !
Jaekwang Kim
--
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.
