This would be for the interpolation only right? So if I use this GlobalRHS, and construct my LHS based on the actual problem I want to solve (just a Poisson, so I would be using dphi), it would be equivalent to solving the problem with the RHS being interpolated?
Karen On Fri, Mar 26, 2010 at 5:23 PM, Rahul Sampath <[email protected]>wrote: > To be more precise: > > GlobalMassMatrix = Sum_over_elements{ > integral_over_current_element_using_quadrature(phi_i*phi_j) } > > GlobalRHS = GlobalMassMatrix*Nodal_F_Vector > > On Fri, Mar 26, 2010 at 5:22 PM, Rahul Sampath <[email protected]> > wrote: > > This is what I meant: > > > > MassMatrix = integral(phi_i*phi_j) > > > > RHS = MassMatrix*Nodal_F_Vector > > > > > > On Fri, Mar 26, 2010 at 5:11 PM, Karen Lee <[email protected]> wrote: > >> Sorry, I think I was probably confused. I guess you just meant that I > can > >> simply use integral of f_i phi_i phi_j as my RHS, and my original LHS as > my > >> LHS, and that would already be effectively an interpolation? Please let > me > >> know if this is the correct way to think about it instead of the long > >> message I sent with lots of code... > >> > >> My 2 questions regarding access to data would still remain: > >> > >> 1) I am still having problems accessing MeshData from the assemble > >> function. I would like to do 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?) > >> > >> 2) 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? > >> > >> Thank you so much!!! > >> Karen > >> > >> > >> On Fri, Mar 26, 2010 at 4:57 PM, Karen Lee <[email protected]> wrote: > >>> > >>> 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
