From: liang cheng <goeasyon@gm...>  20100428 20:09:14

> > > 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 (dim1, 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 } 