John —

I just added TRI3SUBDIVISION to quadrature_trap_2D.C as a case, and I think 
it’s working now. I am still debugging things one step at a time, and one place 
that I got a complaint from is fe_lagrange.C — I may be confusing myself right 
now, but I think it’s complaining because TRI3SUBDIVISION is listed under FIRST 
order along with TRI3, but I think TRI3SUBDIVISION should be FOURTH order. I 
think I am confused because TRI3SUBDIVISION is listed in fe_lagrange.C.

Thank you,
Mike

On Aug 28, 2019, at 12:05 PM, John Peterson 
<jwpeter...@gmail.com<mailto:jwpeter...@gmail.com>> wrote:

Hi Mike,

I think this (using QTrap with TRI3SUBDIVISION elements) should work just fine. 
Did you try it yet? Let us know and we can update the case statements for 
quadrature. There are many places like this in the library where TRI3 and 
TRI3SUBDIVISION elements can be used more or less interchangeably, so there's 
probably a bunch we haven't found yet.

--
John


On Wed, Aug 28, 2019 at 12:33 AM Lee, Jae Ho 
<jaeho...@live.unc.edu<mailto:jaeho...@live.unc.edu>> wrote:
Hello,

I am currently implementing something that makes use of TRI3SUBDIVISION 
elements, and one of the initial things that I want to do is to compute the 
nodal coordinates for each deformed element. I am trying to avoid the use of 
MeshFunction to do this. Of course since subdivision shape functions are not 
interpolatory, I can’t just use the solution of the coordinate system directly 
as the coordinates of the element vertices. So what I am trying to use QTrap 
that uses nodal quadrature points and evaluating the shape functions at the 
quadrature points for a given element to get something like

std::unique_ptr<libMesh::FEBase> fe(libMesh::FEBase::build(2, 
libMesh::SUBDIVISION));
std::unique_ptr<libMesh::QBase> qrule = libMesh::QBase::build(libMesh::QTRAP, 
2, libMesh::FIRST);
fe->attach_quadrature_rule(qrule.get());
fe->reinit(elem);
for (unsigned int qp = 0; qp < qrule->n_points(); ++qp)
{
    for (unsigned int d = 0; d < 3; ++d)
    {
        for (unsigned int i = 0; i < phi.size(); ++i)
        {
            X_node[qp][d] += X_solution[i][d] * phi[i][qp];
        }
    }
}

X_solution is where the solution of the coordinate system has been stored.

Does this seem reasonable?
But quadrature_trap_2D.C complained because TRI3SUBDIVISION is not supported. 
Is it safe to hardcode TRI3SUBDIVISION as a case in the code? Or is there a 
different way to go about it (e.g. is there a different nodal quadrature rule 
that supports TRI3SUBDIVISION)?

Thank you,
Mike



_______________________________________________
Libmesh-users mailing list
Libmesh-users@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/libmesh-users

Reply via email to