Thanks Rahul. Your responses have clarified things for me.
I am still having problems accessing MeshData from the assemble function
(which is only supposed to have 2 arguments right?)
I can do:
void assemble_load(EquationSystems& es,
const std::string& system_name)
{
libmesh_assert (system_name == "load");
const MeshBase& mesh = es.get_mesh();
printf("mesh_data obtained");
......
But I'm not sure how I can get the MeshData object that I created in main
that's attached to the mesh I created. In the look iterating over the
elements el, I know that I have to use something like
mesh_data.get_data(el.get_node(i))[0] to get the first data variable for
node i with an element, i from 0 to 3, but I'm not sure how to access
mesh_data from my System. (I should use Linear Implicit System so I can have
the matrix on the lhs right?)
My code is based on example 3, and the relevant part is this (not sure if
it's correct):
for ( ; el != end_el ; ++el)
{
const Elem* elem = *el;
dof_map.dof_indices (elem, dof_indices);
fe->reinit (elem);
Ke.resize (dof_indices.size(),
dof_indices.size());
Fe.resize (dof_indices.size());
for (unsigned int qp=0; qp<qrule.n_points(); qp++)
{
for (unsigned int i=0; i<phi.size(); i++)
for (unsigned int j=0; j<phi.size(); j++)
{
Ke(i,j) += JxW[qp]*(phi[i][qp]*phi[j][qp]);
}
{
for (unsigned int i=0; i<phi.size(); i++)
{
for (unsigned int k=0; k<phi.size(); k++)
{
Fe(i) +=
JxW[qp]*mesh_data.get_data(el.get_node(k))[0]*phi[i][qp]phi[k][qp];
}
}
}
....
The other question is, for Fe, do I integrate over all the quadrature points
of phi_i and phi_j with f_i being a constant?
Then once I get this solution (the variable name is "R", let's say), I hope
to use it in place of mesh_data.get_data(el.get_node(k))[0] and do
something like this (of course this is another system that I'm adding to
EquationSystems):
*for* (*unsigned* *int* qp=0; qp<qrule.n_points(); qp++)
{
*for* (*unsigned* *int* i=0; i<phi.size(); i++)
*for* (*unsigned* *int* j=0; j<phi.size(); j++)
{
Ke2(i,j) += JxW[qp]*(dphi[i][qp]*dphi[j][qp]);
}
{
*const* Real x = q_point[qp](0);
*const* Real y = q_point[qp](1);
*const* Real z = q_point[qp](2);
*const* Real fxy = R(x,y,z);
*for* (*unsigned* *int* i=0; i<phi.size(); i++)
Fe2(i) += JxW[qp]*fxy*phi[i][qp];
}
}
I'm not sure how to access anything in a Variable object though... Let's say
"R" being a a variable I add to the first system to get the RHS
interpolation, and "u" is the variable I add to the second equation system,
which is the actual solution I'm after... I just know that I can output the
values of the solution at various nodal points in a file, but am not sure
what to do with the data structure and how I can extract values at different
arbitrary locations.
Apologies for the lengthy email...
Thanks,
Karen
On Fri, Mar 26, 2010 at 9:15 AM, Rahul Sampath <[email protected]>wrote:
> Hi Karen:
>
> Take a look at any Nonlinear example problem. Whenever you want to use
> any solution vector in your residual computation, you will need to
> interpolate the nodal values using the FEM shape functions for this
> element and then do the integration. It is very similar to what you
> want to do. That is why I suggested the Mass matrix trick. It is very
> simple to implement and fast too especially if you want to change f
> often. You can use the same Mass matrix with differen f. The catch is
> that you are using your FE shape functions for your interpolation. As
> long as you are happy with linear interpolation for f this should do.
> If you want to interpolate f with a higher order polynomial than your
> FE shape function then this wont work.
>
> Btw, if I was not clear earlier:
> You have to form a global Mass matrix by integrating phi_i phi_j over
> all elements and doing a typically FEM assembly. Then you can simply
> multiply this global Mass matrix with you global nodal vector for f.
>
> On Fri, Mar 26, 2010 at 12:51 AM, Karen Lee <[email protected]> wrote:
> > Hi Rahul, I'm not completely sure what you mean.
> >
> > I would like to form my RHS by integrating f_i Phi_i (I guess there's no
> > need to multiply Phi_j? But you can correct me) for each element.
> >
> > In order to do so, I need values of f at various quadrature points. I
> have f
> > at various nodal values. The question is, how do I get this linear
> > interpolation...
> >
> > Do you mean that, for each element, I form the mass matrix by the xyz
> values
> > of the nodes (and a constant 1) and solve for the coefficient by saying
> > \sum_j A_ij y_j= f_i, where:
> >
> > A = [1 x1 y1 z1;
> > 1 x2 y2 z2;
> > 1 x3 y3 z3;
> > 1 x4 y4 z4] and y_j would be my unknown (where j = 1 corresponds
> to
> > the constant value, and 2, 3, 4 corresponds to the gradient in the x, y,
> z
> > directions respectively)?
> >
> > Thanks,
> > Karen
> >
> >
> > On Thu, Mar 25, 2010 at 11:44 PM, Rahul Sampath <[email protected]
> >
> > wrote:
> >>
> >> If you want to form a RHS by integrating f_i Phi_i Phi_j, You could
> >> form a Mass matrix and then multiply with your vector of nodal values.
> >>
> >> Rahul
> >>
> >> On Thu, Mar 25, 2010 at 11:40 PM, Karen Lee <[email protected]> wrote:
> >> > I'm afraid you misunderstood. I don't have the function that when
> given
> >> > x,
> >> > y, z values gives me the function value. What I do have is just the
> >> > values
> >> > at the nodes of the mesh, which need to be linearly interpolated such
> >> > that I
> >> > will have something like exact_function. which gives me the value when
> >> > supplied with any x, y, z.
> >> >
> >> >
> >> >
> >> > On Thu, Mar 25, 2010 at 10:54 PM, Liang <[email protected]> wrote:
> >> >
> >> >> Karen Lee wrote:
> >> >>
> >> >>> I guess I'm not clear how to do this: Load data as a solution into
> >> >>> that,
> >> >>> and
> >> >>> query
> >> >>> it when you're integrating your real system.
> >> >>>
> >> >>> I have:
> >> >>> Mesh mesh(3);
> >> >>> MeshData mesh_data(mesh);
> >> >>> mesh_data.activate();
> >> >>> mesh.read (mesh_file, &mesh_data);
> >> >>> mesh_data.read(mesh_file);
> >> >>> EquationSystems equation_systems (mesh);
> >> >>>
> >> >>>
> >> >>> equation_systems.add_system<ExplicitSystem> ("RHS");
> >> >>> equation_systems.get_system("RHS").add_variable("R", FIRST);
> >> >>>
> >> >>> After that, I'm not clear how exactly to load data as a solution in
> >> >>> the
> >> >>> code. My goal is to get a linearly interpolated function of my data
> on
> >> >>> the
> >> >>> nodes (in the form of exact_solution, such that I get the function
> >> >>> value
> >> >>> out
> >> >>> when supplying x, y and z).
> >> >>>
> >> >>> Hope that clarifies things, and sorry for the multiple emails...
> >> >>>
> >> >>> Karen
> >> >>>
> >> >>>
> >> >>>
> ------------------------------------------------------------------------------
> >> >>> Download Intel® Parallel Studio Eval
> >> >>> Try the new software tools for yourself. Speed compiling, find bugs
> >> >>> proactively, and fine-tune applications for parallel performance.
> >> >>> See why Intel Parallel Studio got high marks during beta.
> >> >>> http://p.sf.net/sfu/intel-sw-dev
> >> >>> _______________________________________________
> >> >>> Libmesh-users mailing list
> >> >>> [email protected]
> >> >>> https://lists.sourceforge.net/lists/listinfo/libmesh-users
> >> >>>
> >> >>>
> >> >>>
> >> >> so you already have the function, which is obtained from your
> discreted
> >> >> data?
> >> >> then just put the function as the exact_function.
> >> >> I think you are trying the 3D case, start from a 2d will be easier.
> >> >>
> >> >> Liang
> >> >>
> >> >
> >> >
> ------------------------------------------------------------------------------
> >> > Download Intel® Parallel Studio Eval
> >> > Try the new software tools for yourself. Speed compiling, find bugs
> >> > proactively, and fine-tune applications for parallel performance.
> >> > See why Intel Parallel Studio got high marks during beta.
> >> > http://p.sf.net/sfu/intel-sw-dev
> >> > _______________________________________________
> >> > Libmesh-users mailing list
> >> > [email protected]
> >> > https://lists.sourceforge.net/lists/listinfo/libmesh-users
> >> >
> >
> >
>
------------------------------------------------------------------------------
Download Intel® Parallel Studio Eval
Try the new software tools for yourself. Speed compiling, find bugs
proactively, and fine-tune applications for parallel performance.
See why Intel Parallel Studio got high marks during beta.
http://p.sf.net/sfu/intel-sw-dev
_______________________________________________
Libmesh-users mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/libmesh-users
