Creating a new Element (some guide lines, please)

[Takala]

Hello,

I intend to create a new element to Elmer.

In my case I need to add equations and variables to the Whitney element. In addition to A and V, I need to have induction components b2, b4, etc. for modeling a lamination stack with a homogenization technique.

Please give some guidelines how to do this. For example what are the files, functions and subroutines to look at when creating an element (or in this case when modifying the Whitney element)?

Best Regards,

Eelis

[Takala]

So I have been digging around and now I have some hunch think where to go:

I will make a test case where the solution variable for the element should be the following

(v0, v1, ... v_np, A0, A1, ..., A_nd, b2)

where the b2 is the new variable which is describing the average induction component (parallel to the lamination) averaged over the element.

now, I have seen that the "Variable % Dofs" is 1 currently for the MagnetoDynamics solver and the element type is "n:1 e:1".
I think that this means that the Variable is one dimensional and thus has one degree of freedom for node and edge. For example
for tetraeder element the solution vector is (v0, v1, v2, v3, A0, A1, A2, A3, A4, A5, A6). Now I need to find a way to tell to Elmer
that there should be one more degree of freedom for the b2, for example the tetraeder would have (v0, v1, v2, v3, A0, A1, A2, A3, A4, A5, A6, b2)

How could this be done?

Should I come up with some Element code "n:1 e:1 ind:1" and somehow change the element accordingly. If yes, how to do this?

Best Regards,

Eelis

[Takala]

OK, I get it now.

I could use the BDOFS to increase the number of degrees of freedom in the element:

"MainUtils.src":

Code: Select all

MaxBDOFs = 0
          DO i=1,Solver % Mesh % NumberOFBulkElements
            CurrentElement => Solver % Mesh % Elements(i)
            MaxBDOFs  = MAX( MaxBDOFs,  CurrentElement % BDOFs )
          END DO

This could be done by saying

Code: Select all

Element = "n:1 e:1 b:1"

Right?

One question remains:
when using the Element keyword in the solver section, the keyword is applied to all the elements activated for the specific solver. What if one would like to apply to only a subset of the active elements? Clearly the element type structure and the variable initialization routine supports this. But is there an easy way to do this? Some sif keyword?

Cheers,

Eelis

[raback]

Hi Eelis

So you would like to activate some dofs in the element definition only for some bodies? I don't think there is an elegant existing way of doing this. There is a branch for doing this to nodal elements but not for these generic elements. The ugly remedy is to set in the non-active bodies for given dofs some dummy equation in the style that Dirichlet conditions are set.

Could you ellaborate a little on your equation. What are these induction components b2, b4 etc.?

-Peter

[Takala]

Hi Peter,

yes that's right, I would like to add dofs to elements in order to describe a lamination stack without drawing the laminates. Would it be a lot of work to have the dofs only for some elements? How does it work for nodal elements? What spell you use in sif for that?

Well, even though there isn't a good way to add the dofs for specified bodies in the case of generic elements, it doesn't stop me from testing the equations. It just slows things down I guess. But now there is a demand, so there is a reason to provide, right? Just kidding, I know you are busy.

The homogenization technique that I'm following is described in [1]. The point is to describe the behavior of a steel laminate stack without meshing every single laminate. The conventional weak form of the Ampere's law is
(v curl a, curl a') + (sigma d_t(a), a') = (j,a')

If the so called low field approximation is used, one can just add a term (sigma*d_l^2/12 d_t(curl a), curl a') to the equation. However, this approximation is only for very low frequencies. Now in [1] the so called second order approximation is used, there is a new dof called b2 which is a component that describes the average flux in plane of the steel laminate. Similarly if a higher order is chosen, the dofs (induction components, b3, b4..) increase. These components are coming from the polynomial expansion of the magnetic flux density. There is another expansion in Paavo Rasilo's thesis where Fourier series are used. But essentially it is the same technique. Thus when these dofs are introduced, there are equations introduced accordingly, for example the second order approximation has two equations:

(v curl, a, curl a') +(v_t curl a, curl a') + (sigma_{perp} d_t(a), a') + ( (sigma_{par} d_l^2)/12 * d_t(curl a), curl a') - ( (sigma_{par} d_l^2)/12 * d_t(curl a), curl a') = (j,a')

and

(v_t/5*b2, b2') - ( (sigma_{par} d_l^2)/60 * d_t(curl a), b2')+(sigma_{par} d_l^2/210 d_t(b2), b2')

where v_t, sigma_{perp} and sigma_{par} are the reluctivity, perpendicular conductivity and parallel conductivity tensors, which have 3, 2 and 2 non-zero components, respectively.
when the third order is used there are three equations, etc.

By the way, about the tensors: the technical implementation of getting tensors from sif to the f90 module is done for example in "HeatSolve.src" in this way:

Code: Select all

CALL ListGetRealArray( Material,'Heat Conductivity',Hwrk,n, &
                      Element % NodeIndexes )
         HeatConductivity = 0.0d0
         IF ( SIZE(Hwrk,1) == 1 ) THEN
           DO i=1,3
             HeatConductivity( i,i,1:n ) = Hwrk( 1,1,1:n )
           END DO
         ELSE IF ( SIZE(Hwrk,2) == 1 ) THEN
           DO i=1,MIN(3,SIZE(Hwrk,1))
             HeatConductivity(i,i,1:n) = Hwrk(i,1,1:n)
           END DO
         ELSE
           DO i=1,MIN(3,SIZE(Hwrk,1))
             DO j=1,MIN(3,SIZE(Hwrk,2))
               HeatConductivity( i,j,1:n ) = Hwrk(i,j,1:n)
             END DO
           END DO
         END IF

Is this the way to do it? Just a curiosity: where does the Hwrk name come from?

[1] J. Gyselink and P. Dular, "A Time-Domain Homogenization Technique for Laminated Iron Cores in 3-D Finite-Element Models", IEEE Trans. Magn. 40(2) March 2004

[mika]

Hi,

I think there is a trick that enables using bodywise element definitions even when the element type definition is not classic. The following worked for me a couple of years ago with hp-approximation (and I hope still does), so you might consider trying it.

For example, if we have split the entire computational domain into two bodies and we want to use basis functions of degree p = 1 in the body 1 and the polynomial degree p = 3 in the body 2, we can make the desired element choices by first making two copies of the equation section. These two equation sections thus contain the same solvers list. We may then differentiate them by giving different element definitions. The whole example construction in the case of the solvers list of the size 2 is as follows:

Body 1
Equation = 1
Material = 1
End

Body 2
Equation = 2
Material = 1
End

Equation 1
Active Solvers(2) = 1 2
Element{1} = String "p:1"
Element{2} = String "p:1"
End

Equation 2
Active Solvers(2) = 1 2
Element{1} = String "p:3"
Element{2} = String "p:3"
End

Note that the integer identifier between the braces {} refers to the integer identifier of the solver section. In your case, you might consider replacing the element type strings "p:1" and "p:3" by "n:1 e:1" and "n:1 e:1 b:1", for example. When I was working with this, I tested this construction only in the case of p-elements. Whether this works when the number of edge (e) and bubble/elementwise (b) degrees of freedom is additionally specified may have not been tested in a thorough manner or even at all.

Best regards,
Mika

[Takala]

Hi Mika,

very nice. I think this might actually be perfect for me if it does work.
I will report here the feasibility of this trick as soon as I have time to test it. I still need to do the equations first.

btw. I tested the low field approximation and it works well and is easy to make.

Cheers,

Eelis