>
>
> And what about the boundary fe object? Did you also reinit that and
> with the correct side?
>
> --
> John
>
Yes, I did the BC like the ex8, the Neumann BC is defined in the assemble()
and
the Dirichlet BC is defined in fill_dirichlet_bc().
Is it critical if no fe_face->attach_quadrature_rule() and reinit() in the
fill_dirichlet_bc()?
No quadrature rule was added in the fill_dirichlet_bc() because no such
commands in the
ex8 neither, I don't dare to violate the example too much, so didn't plug
the quadrature
rule in fill_dirichlet_bc().
I guess this problem is related to the matrix definition, maybe my guessing
is incorrect.
The matrix system was defined as followed,
SparseMatrix<Number>& stiffness =
large_deform_system.get_matrix("stiffness");
SparseMatrix<Number>& damping =
large_deform_system.get_matrix("damping");
SparseMatrix<Number>& mass = large_deform_system.get_matrix("mass");
NumericVector<Number>& force =
large_deform_system.get_vector("force");
However these matrix doesn't participate in computation. These matrix seems
required for the Newmark
solver as I dig into the newmark_system.C/h. Is this 'SparseMatrix' data
structure
quit different with the that of 'DenseMatrix<Number>'? I found the
'DenseSubMatrix'
could be defined while the 'SparseSubMatrix' not. Since the governing
equation here
require the matrix structure like
DenseSubMatrix<Number>
Kuu(Ke), Kuv(Ke), Kuw(Ke),
Kvu(Ke), Kvv(Ke), Kvw(Ke),
Kwu(Ke), Kwv(Ke), Kww(Ke);
So the question is building Sub Matrix for the SparseMatrix like 'stiffness
'.
Maybe I mix the two separated problems. I have double checked the codes
with the ex8, the only difference up to now is the submatrix system. For the
variable of wave equation in ex8 is only one (pressure), while here in my
case
is three variables, which leads to a 3x3 stiffness matrix.
I thought the solid element is easy, but it turns out almost kill me. I
appreciate
your time and answer. Thanks!
Liang
********************************************************************
{ // boundary condition below
for (unsigned int s=0; s<elem->n_sides(); s++)
if (elem->neighbor(s) == NULL) {
* AutoPtr<FEBase> fe_face (FEBase::build(dim, fe_type));
QGauss qface (dim-1, fe_type.default_quadrature_order());
fe_face-> attach_quadrature_rule (&qface);*
const std::vector<std::vector<Real> >& phi_face =
fe_face->get_phi();
const std::vector<Real>& JxW_face = fe_face->get_JxW();
* fe_face->reinit(elem,s);*
for (unsigned int qp=0; qp<qface.n_points(); qp++)
{
// Neumman boundary condition
for (unsigned int i=0; i<phi_face.size(); i++){
Fu(i) += JxW_face[qp] * q_n * phi_face[i][qp];
Fv(i) += JxW_face[qp] * q_n * phi_face[i][qp];
Fw(i) += JxW_face[qp] * q_n * phi_face[i][qp];
}
} // end face node loop
} // end if (elem->neighbor(side) == NULL)
} // end BC
/*******************************************************************
********************************************************************
*** Apply the dirichlet boundary condition u, v , w ***
********************************************************************
*******************************************************************/
void fill_dirichlet_bc(EquationSystems& es,
const std::string& system_name)
{
const Real y_coo = 6.; //BC coordinates
libmesh_assert (system_name == "Disp_finite_deform_3d");
NewmarkSystem & system = es.get_system<NewmarkSystem> (system_name);
SparseMatrix<Number>& matrix = *system.matrix;
NumericVector<Number>& rhs = *system.rhs;
const MeshBase& mesh = es.get_mesh();
const bool do_for_matrix =
es.parameters.get<bool>("Newmark set BC for Matrix");
const Real current_time = es.parameters.get<Real>("time");
unsigned int n_nodes = mesh.n_nodes();
for (unsigned int n_cnt=0; n_cnt<n_nodes; n_cnt++)
{
const Node& curr_node = mesh.node(n_cnt);
const Real penalty = 1.e10; // The penalty parameter.
if (fabs(curr_node(2) - 0.0) < TOLERANCE)
{
unsigned int dn = curr_node.dof_number(0,0,0);
Real u_value = .0;
Real v_value = .0;
Real w_value = .0;
rhs.add(dn, u_value*penalty);
rhs.add(dn, v_value*penalty);
rhs.add(dn, w_value*penalty);
if (do_for_matrix)
matrix.add(dn, dn, penalty);
}
// force or displacement
if (fabs(curr_node(2)-y_coo) < TOLERANCE)
{
unsigned int dn = curr_node.dof_number(0,0,0);
Real w_value;
if (current_time < .002 )
w_value = -0.1; // sin(2*pi*current_time/.002);
else
w_value = -0.2;
rhs.add(dn, w_value*penalty);
if (do_for_matrix)
matrix.add(dn, dn, penalty);
}
} // loop n_cnt
}
------------------------------------------------------------------------------
_______________________________________________
Libmesh-users mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/libmesh-users
