## Re: [Libmesh-users] XYZ spatial coordinates

 Re: [Libmesh-users] XYZ spatial coordinates From: Roy Stogner - 2008-03-04 14:16:05 ```On Tue, 4 Mar 2008, Maxime Debon wrote: > Looking for tips to check the correctness of Ke matrix and Fe vector > in an handbook, I found this rule : A0 * u'' = (A0 * u')' - A0' * u' > > This latest transformation must be taken into account to adjust a > non-constant coefficient to the diffusion term. It also explains why > the result was correct with a constant coefficient. Oh! I didn't even think of that. Yeah; if the coefficient is outside both derivatives then you get an extra term from it when integrating by parts. I'm more used to terms like (A0*u')' (which doesn't give you an extra term) than I am to A0*u''. Double-check and make sure that the latter is really what models your physical system. Often conservation models have a conservation statement: D/dt(conserved) = div(flux) and a constitutive relation: flux = A*grad(conserved) In which case A really does end up inside one but not both derivatives in the strong form of the resulting equation. --- Roy ```

 [Libmesh-users] Newmark System - 1D From: Maxime Debon - 2008-02-17 11:00:43 ```Hi, As the 2D heat equation's results are quite good under the newmark system, I tried to implement the same equation in 1D to compare with an analytical solution. Not convinced by the solution values in 1D, I found that the initial condition (set to 1) was still 0 after assembling the system. Furthermore, even if the project solution tool gave satisfaction to affect initial values with 1, the final result kept going ambiguous. I am almost sure to make a mistake but, does the newmark class is configured to work with 1D meshes ? Best regards, Maxime NB: Libmesh is configured with PETSc and the code comes from example 8 ```
 [Libmesh-users] XYZ spatial coordinates From: Maxime Debon - 2008-02-27 19:48:16 ```Hi, Using spatial coordinates (q_point) in matrix Ke almost as it's presented in the example 14 for matrix Fe, I may misunderstand something important. Here is the code. // ------------------------------------------------------------------ // const std::vector& qrule_points = fe->get_xyz(); for (unsigned int qp=0; qp
 Re: [Libmesh-users] XYZ spatial coordinates From: Maxime Debon - 2008-03-03 22:44:27 ```Hi, To give a precision, the variable coefficient only generate an error computing the diffusion term (u''(x)). In contrast, tests made with an X dependent coordinate coefficient on convection term and u(x) term gave a very powerful precision. Thus, the error can't be attributed to previous definitions in the code. // ----------------------------------------------------------------------- // Taking the equation : x * u''(x) = 1 between [0.5 ; 1.0] with u[0.5] = 0.5 and u[1.0] = 1.0, we would expect a value of 0.88045 at the point 0.900391 but we obtain 0.900391. The corresponding code is : Ke(i,j) += JxW[qp]*(dphi[i][qp]*dphi[j][qp])*(-q_point[qp](0)); Fe(i) += JxW[qp]*phi[i][qp]*1.0; The correct result could be obtained computing " u''(x) = 1/x ", making disappear the variable coefficient on the diffusion term. The corresponding code is : Ke(i,j) += JxW[qp]*(dphi[i][qp]*dphi[j][qp])*(-1.0); Fe(i) += JxW[qp]*phi[i][qp]*1.0*q_point[qp](0); However this kind of solution is only available in an 1D equation. // ----------------------------------------------------------------------- // These results are obtained with the linear system presented in example 0. Furthermore, the problem only occurs with a non-constant diffusion term coefficient. Ideas and comments welcome, Best regards, Maxime ```
 Re: [Libmesh-users] XYZ spatial coordinates From: Roy Stogner - 2008-03-03 22:56:42 ```On Mon, 3 Mar 2008, Maxime Debon wrote: > Taking the equation : x * u''(x) = 1 > between [0.5 ; 1.0] > with u[0.5] = 0.5 and u[1.0] = 1.0, > > we would expect a value of 0.88045 at the point 0.900391 but we obtain > 0.900391. This problem is still giving you trouble? I recall seeing your earlier post but figured it was probably just a bug in user code that you'd have found yourself before I could send an email. Some questions: Is the bad solution you're getting just a straight line from 0.5,0.5 to 1,1? Or is the value 0.900391,0.900391 just a coincidence? What kind of mesh, finite element type, polynomial degree, etc are you using? Have you tried verifying that the element matrices and load vectors are correct? For a 1D problem with p=1 this is very easy. > The corresponding code is : > Ke(i,j) += JxW[qp]*(dphi[i][qp]*dphi[j][qp])*(-q_point[qp](0)); > Fe(i) += JxW[qp]*phi[i][qp]*1.0; > > > The correct result could be obtained computing " u''(x) = 1/x ", > making disappear the variable coefficient on the diffusion term. > > The corresponding code is : > Ke(i,j) += JxW[qp]*(dphi[i][qp]*dphi[j][qp])*(-1.0); > Fe(i) += JxW[qp]*phi[i][qp]*1.0*q_point[qp](0); Oddly enough; I'd have expected just the opposite behavior. If you're working from example 0, the 5th order quadrature rule along with 2nd order elements should integrate u'*v'*x perfectly, but should give you some quadrature error with 1/x*v. --- Roy ```
 Re: [Libmesh-users] XYZ spatial coordinates From: Maxime Debon - 2008-03-04 13:49:36 ```Hi, Looking for tips to check the correctness of Ke matrix and Fe vector in an handbook, I found this rule : A0 * u'' = (A0 * u')' - A0' * u' This latest transformation must be taken into account to adjust a non-constant coefficient to the diffusion term. It also explains why the result was correct with a constant coefficient. Thank you for your advice, Best regards, Maxime ```
 Re: [Libmesh-users] XYZ spatial coordinates From: Roy Stogner - 2008-03-04 14:16:05 ```On Tue, 4 Mar 2008, Maxime Debon wrote: > Looking for tips to check the correctness of Ke matrix and Fe vector > in an handbook, I found this rule : A0 * u'' = (A0 * u')' - A0' * u' > > This latest transformation must be taken into account to adjust a > non-constant coefficient to the diffusion term. It also explains why > the result was correct with a constant coefficient. Oh! I didn't even think of that. Yeah; if the coefficient is outside both derivatives then you get an extra term from it when integrating by parts. I'm more used to terms like (A0*u')' (which doesn't give you an extra term) than I am to A0*u''. Double-check and make sure that the latter is really what models your physical system. Often conservation models have a conservation statement: D/dt(conserved) = div(flux) and a constitutive relation: flux = A*grad(conserved) In which case A really does end up inside one but not both derivatives in the strong form of the resulting equation. --- Roy ```