| Dear Johnathan,
The expressions I gave to you are correct, you can enforce the symmetry M=gf_asm('volumic','M(#1,#1)+=sym(comp(vBase(#1).vBase(#1))(:,k,:,k))', mim,mf);
K=gf_asm('volumic','M(#1,#1)+=sym(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);
Then you need to calculate the eigen vectors of M^-1 K Be aware that the Kernel is very big. I cite Ronan: "Your matrix K should be positive semi-definite with a large kernel equal to the number of nodes in your mesh minus one (modulo the boundary conditions)."
Best Regards
Louis
Dear Johnathan,
When using finite element, you are solving the governing equation (in your case, the Maxwell equation) using a weak formulation. The discretisation of the unknown field injected into the weak formulation leads to an algebraical system. It seems to me that you're not using a weak formulation. The weak formulation for maxwell problem is: " find E in H(curl) such as for any E' in H(curl) \int_\Omega -omega/c^2 E.E' + curlE curlE' dx=0 "
the first part correspond to the "mass matrix" M, the second to the "stiffness matrix" K you get the mode by taking the eigen vector of (M^-1 K).
I have the following code, but unfortunatelly the stiffness matrix doesn't work. I am still trying to find out why.
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)')); M=gf_asm('volumic','M(#1,#1)+=comp(vBase(#1).vBase(#1))(:,k,:,k)', mim,mf);
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);
I hope it will help you.
Best Regards
Louis Kovalevsky
Hello,
I have been trying to calculate the eigenmodes of light travelling through a dielectric waveguide. To do this, I have been using the assembly tools for arbitrary matrices. However, when I attempt to solve for eigenvectors of the stiffness matrix using ARPACK, I do not get sensible results.
I have been attempting to generate matrices for the Maxwell Wave equation using:
Del^2 Psi - (omega^2 / c^2 ) Psi
For the first term, I simply use the laplacian given in the getfem documentation: M(#1, #1)+=comp(Grad(#1).Grad(#1))(:,k,:,k). However, I am unsure as to how I should generate the second term. However, I am not sure how to generate the second term in the equation. Previously, I have used: (comp(Base(#1).Base(#1).Base(#2))(:,:,j).h(j)). Another technique I have considered was to take the vector generated by interpolating the data set representing (omega^2 / c^2 ), converting the elements of the vector to elements of a diagonal matrix, and simply adding the result to the matrix of the first term. Could anyone tell me if either of these approaches is correct, or if the correct procedure is some method that I have not yet thought of?
Yours,
Johnathan Yik
Getfem-users mailing list
address@hidden
https://mail.gna.org/listinfo/getfem-users
Getfem-users mailing list
address@hidden
https://mail.gna.org/listinfo/getfem-users
|