I am working on solving the level set equation using DG methods, form the
looks of example Miscellaneous Example 5 this should be possible.
Question 1:
Currently, I am simply trying to formulate the mass matrix, which is simply
the integral of N^T N over the domain, where N is the shape function
vector. I would like this matrix to be diagonal so it can be inverted
locally allowing me to solve the system with an explicit scheme. Does
libMesh have a set of basis functions that could offer me his behavior? The
book I am following for this problem recommends "orthogonal hierarchical
shape functions". So, I tried to using L2_HIERARCHIC via the
add_system(...) command, but the mass matrix was not diagonal. Could it be
a problem with how I am assembling the matrix via quadrature points, my
assemble() function is shown at the end. I am new to FEM and mostly
self-taught so I am naive in many aspects.
Question 2:
Is there any professional help available for developing libMesh code. I
have a research grant and could offer payment for assistance, which I need
and will need as I continue on this project. There is also the potential
for collaboration and publications if anyone is interested.
Thanks and I appreciate any help I can get.
Andrew
---- ASSEMBLE FUNCTION
------------------------------------------------------------------------
void LevelSetSystem :: assemble(){
// Get a constant reference to the mesh object.
const MeshBase& mesh = get_mesh();
// The dimension that we are running
const unsigned int dim = mesh.mesh_dimension();
// Get a constant reference to the Finite Element type
// for the first (and only) variable in the system.
FEType fe_type = variable_type(0);
// Build a Finite Element object of the specified type.
AutoPtr<FEBase> fe (FEBase::build(dim, fe_type));
AutoPtr<FEBase> fe_face (FEBase::build(dim, fe_type));
// A Gauss quadrature rule for numerical integration.
// Let the \p FEType object decide what order rule is appropriate.
QGauss qrule (dim, fe_type.default_quadrature_order());
QGauss qface (dim-1, fe_type.default_quadrature_order());
// Tell the finite element object to use our quadrature rule.
fe->attach_quadrature_rule (&qrule);
fe_face->attach_quadrature_rule (&qface);
// Here we define some references to cell-specific data that
// will be used to assemble the linear system. We will start
// with the element Jacobian * quadrature weight at each integration point.
const std::vector<Real>& JxW = fe->get_JxW();
const std::vector<Real>& JxW_face = fe_face->get_JxW();
// The element shape functions evaluated at the quadrature points.
const std::vector<std::vector<Real> >& N = fe->get_phi();
const std::vector<std::vector<Real> >& N_face = fe_face->get_phi();
// The element shape function gradients evaluated at the quadrature
// points.
const std::vector<std::vector<RealGradient> >& B = fe->get_dphi();
// The XY locations of the quadrature points used for face integration
const std::vector<Point>& qface_points = fe_face->get_xyz();
// A reference to the \p DofMap object for this system.
const DofMap& dof_map = get_dof_map();
// Define data structures to contain the element matrix
// and right-hand-side vector contribution.
DenseMatrix<Number> Me; // element mass matrix
//DenseMatrix<Number> Ke; // element stiffness matrix
// Storage for the degree of freedom indices
std::vector<unsigned int> dof_indices;
// Get the system time
Real t = this->time;
// Loop over all the elements in the mesh that are on local processor
MeshBase::const_element_iterator el = mesh.active_local_elements_begin();
const MeshBase::const_element_iterator end_el =
mesh.active_local_elements_end();
for ( ; el != end_el; ++el){
// Pointer to the element current element
const Elem* elem = *el;
// Get the degree of freedom indices for the current element
dof_map.dof_indices(elem, dof_indices);
// Compute the element-specific data for the current element
fe->reinit (elem);
// Extract a vector of quadrature x,y,z coordinates
const vector<Point> qp_vec = fe->get_xyz();
// Zero the element matrices and vectors
Me.resize (dof_indices.size(), dof_indices.size());
//Ke.resize (dof_indices.size(), dof_indices.size());
// Compute the RHS and mass and stiffness matrix for this element (Me)
for (unsigned int qp = 0; qp < qrule.n_points(); qp++){
// Extract the velocity vector for the current point
VectorValue<Number> v(dim);
velocity(v, qp_vec[qp], t);
for (unsigned int i = 0; i < N.size(); i++){
for (unsigned int j = 0; j < N.size(); j++){
Me(i,j) += JxW[qp]*(N[i][qp] * N[j][qp]); // mass matrix
//Ke(i,j) += JxW[qp]*(N[i][qp] * (v*B[j][qp])); // stiffness matrix
}
}
}
printf("Me:\n");
Me.print(std::cout);
// Apply the local components to the global mass matrix
request_matrix("mass")->add_matrix(Me, dof_indices);
//this->matrix->add_matrix(Ke, dof_indices);
} // (end) for ( ; el != end_el; ++el)
// Indicate that the mass matrix is complete
request_matrix("mass")->close();
//this->matrix->close();
} // (end) assemble()
--
Andrew E. Slaughter, PhD
[email protected]
Materials Process Design and Control Laboratory
Sibley School of Mechanical and Aerospace Engineering
169 Frank H. T. Rhodes Hall
Cornell University
Ithaca, NY 14853-3801
(607) 229-1829
http://aeslaughter.github.com/
------------------------------------------------------------------------------
Live Security Virtual Conference
Exclusive live event will cover all the ways today's security and
threat landscape has changed and how IT managers can respond. Discussions
will include endpoint security, mobile security and the latest in malware
threats. http://www.accelacomm.com/jaw/sfrnl04242012/114/50122263/
_______________________________________________
Libmesh-users mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/libmesh-users
