|
From: | Ioannis Koufogiannis |
Subject: | Re: [Getfem-users] use of Nedelec element in 3D electromagnetism |
Date: | Wed, 27 Nov 2013 16:28:23 +0100 |
User-agent: | Mozilla/5.0 (Windows NT 6.1; WOW64; rv:24.0) Gecko/20100101 Thunderbird/24.1.0 |
Dear Louis,
some time ago I developed such a code and if I remember well it was working properly. The piece of code that I used for assembling the stiffness matrix, using the C++ getfem functions, was: getfem::generic_assembly matrix_A; ... ... std::stringstream sstr; sstr << "t=comp(vGrad(#1).vGrad(#1));""M(#1,#1)+=sym((t(:,j,k,:,j,k) - t(:,j,k,:,k,j) - t(:,k,j,:,j,k) + t(:,k,j,:,k,j))*0.5)"; matrix_A.set(sstr.str()); I arrived to this result after some expansions and rearrangements of the vectorial products on the paper. At the moment, I am not sure how this differs from your development, but you can try it. It seems that in my calculations, compared to yours, I have all the combinations of (j,k), 1<j,k<3 and then I half the final result. Best, Yiannis On 27.11.2013 14:27, Louis Kovalevsky wrote: Dear Getfem++ users, I would like to use Nedelec finite element to calculate the mode of 3d electromagnetic cavity. I'm using the matlab interface. The importation of the mesh, definition of element and integration method is fine: m = gf_mesh('import', 'gmsh', 'cylinder.msh'); mf = gf_mesh_fem(m,3); gf_mesh_fem_set(mf,'fem',gf_fem('FEM_NEDELEC(3)')); mim=gfMeshIm(m,gf_integ('IM_TETRAHEDRON(8)')); To assemble the term \int E.E'dx I'm using: Ma=gf_asm('volumic','M(#1,#1)+=comp(vBase(#1).vBase(#1))(:,k,:,k)', mim,mf); I got the same result with "gf_asm('mass matrix', mim, mf);" ( as excepted) For the "stiffness" matrix \int curlE curlE' dx i'm using: K=gf_asm('volumic','M(#1,#1)+=comp(vGrad(#1).vGrad(#1))(:,2,3,:,2,3)+comp(vGrad(#1).vGrad(#1))(:,3,2,:,3,2)-comp(vGrad(#1).vGrad(#1))(:,2,3,:,3,2)-comp(vGrad(#1).vGrad(#1))(:,3,2,:,2,3)+comp(vGrad(#1).vGrad(#1))(:,3,1,:,3,1)+comp(vGrad(#1).vGrad(#1))(:,1,3,:,1,3)-comp(vGrad(#1).vGrad(#1))(:,3,1,:,1,3)-comp(vGrad(#1).vGrad(#1))(:,1,3,:,3,1)+comp(vGrad(#1).vGrad(#1))(:,1,2,:,1,2)+comp(vGrad(#1).vGrad(#1))(:,2,1,:,2,1)-comp(vGrad(#1).vGrad(#1))(:,1,2,:,2,1)-comp(vGrad(#1).vGrad(#1))(:,2,1,:,1,2)', mim,mf); My problem is that the obtained matrix K is not positive definite (and it should be...). Does anyone have a suggestion why the K matrix is not positive definite? Thank you in advance for the answer. Best Regards Louis Kovalevsky |
[Prev in Thread] | Current Thread | [Next in Thread] |