If you can tell me how your f_value is stored in your files, I can tell you how to form F_nodal.
On Fri, Mar 26, 2010 at 5:53 PM, Rahul Sampath <[email protected]> wrote: > I have not used MeshData and may be some of the LibMesh experts in the > group can help you with that. However, all you need is a Petsc Vec > object for F_nodal since you can simply call Petsc MatMult to form RHS > once you have Petsc Mat for MassMatrix and F_nodal. > If you know the global DOF id for each mesh node and the corresponding > f_value then you can set values into F_nodal using Petsc functions. > > On Fri, Mar 26, 2010 at 5:43 PM, Karen Lee <[email protected]> wrote: >> Yeah, this makes perfect sense to me. >> >> I'm having difficulty with step 2 though. In main, I define a MeshData that >> I read in from an input file with the Mesh itself and the associated data. >> Basically, F_nodal IS mesh_data.get_data(el.get_node(i))[0], but I don't see >> a get_mesh_data function or something similar in the class doc for System... >> Nor can I find a way to access the MeshData object associated with the mesh >> once I have the mesh object pointer... >> >> >> >> >> On Fri, Mar 26, 2010 at 5:38 PM, Rahul Sampath <[email protected]> >> wrote: >>> >>> Here are the steps to solve Laplacian u(x,y,z) = f(x,y,z): >>> >>> 1. Form A = Assembly( A_elem ) where A_elem = integral(dPhi_i*dPhi_j) >>> 2. Form F_nodal = f(x_i,y_i,z_i) at all mesh nodes >>> 3. Form MassMatrix = Assembly( M_elem ) where M_elem = integral(Phi_k, >>> Phi_j) >>> 4. Form RHS = MassMatrix*F_nodal >>> 5. Solve A U_nodal = RHS for U_nodal >>> >>> >>> >>> On Fri, Mar 26, 2010 at 5:26 PM, Karen Lee <[email protected]> wrote: >>> > 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
