bele.c

Go to the documentation of this file.

/*
    This file is part of magpar.

    Copyright (C) 2002-2009 Werner Scholz

    www:   http://www.magpar.net/
    email: magpar(at)magpar.net

    magpar is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 2 of the License, or
    (at your option) any later version.

    magpar is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with magpar; if not, write to the Free Software
    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
*/

#include "field.h"
#include "util/util.h"

/*
  hybrid FEM/BEM method from:
  T. R. Koehler, D. R. Fredkin, "Finite Element Methods for
  Micromagnetics", IEEE Trans. Magn. 28 (1992) 1239-1244

  matrix elements calculated according to:
  D. A. Lindholm, "Three-Dimensional Magnetostatic Fields from
  Point-Matched Integral Equations with Linearly Varying Scalar
  Sources", IEEE Trans. Magn. MAG-20 (1984) 2025-2032.
*/

int Bele(PetscReal* bvert,PetscReal* facv1,PetscReal* facv2,PetscReal* facv3,PetscReal* matele)
{
  PetscReal   *rr,zeta[ND],zetal;
  PetscReal   rho1[ND],rho2[ND],rho3[ND];
  PetscReal   rho1l,rho2l,rho3l;
  PetscReal   s1,s2,s3;
  PetscReal   eta1[ND],eta2[ND],eta3[ND];
  PetscReal   eta1l,eta2l,eta3l;
  PetscReal   xi1[ND],xi2[ND],xi3[ND];
  PetscReal   gamma1[ND],gamma2[ND],gamma3[ND];
  PetscReal   p[ND],a,omega;

  matele[0]=matele[1]=matele[2]=0.0;

  /* get coordinates of face's vertices */
  my_dcopy(ND,facv1,1,rho1,1);
  my_dcopy(ND,facv2,1,rho2,1);
  my_dcopy(ND,facv3,1,rho3,1);

  /* calculate edge vectors and store them in xi_j */
  my_dcopy(ND,rho2,1,xi1,1);
  my_daxpy(ND,-1.0,rho1,1,xi1,1);
  my_dcopy(ND,rho3,1,xi2,1);
  my_daxpy(ND,-1.0,rho2,1,xi2,1);
  my_dcopy(ND,rho1,1,xi3,1);
  my_daxpy(ND,-1.0,rho3,1,xi3,1);

  /* calculate zeta direction */
  douter(ND,xi1,xi2,zeta);

  /* calculate area of the triangle */
  zetal=my_dnrm2(ND,zeta,1);
  a=0.5*zetal;

  /* renorm zeta */
  my_dscal(ND,1.0/zetal,zeta,1);

  /* calculate s_j and normalize xi_j */
  s1=my_dnrm2(ND,xi1,1);
  my_dscal(ND,1.0/s1,xi1,1);
  s2=my_dnrm2(ND,xi2,1);
  my_dscal(ND,1.0/s2,xi2,1);
  s3=my_dnrm2(ND,xi3,1);
  my_dscal(ND,1.0/s3,xi3,1);

  douter(ND,zeta,xi1,eta1);
  douter(ND,zeta,xi2,eta2);
  douter(ND,zeta,xi3,eta3);

  gamma1[0]=gamma3[1]=my_ddot(ND,xi2,1,xi1,1);
  gamma1[1]=my_ddot(ND,xi2,1,xi2,1);
  gamma1[2]=gamma2[1]=my_ddot(ND,xi2,1,xi3,1);

  gamma2[0]=gamma3[2]=my_ddot(ND,xi3,1,xi1,1);
  gamma2[2]=my_ddot(ND,xi3,1,xi3,1);

  gamma3[0]=my_ddot(ND,xi1,1,xi1,1);

  /* get R=rr */
  rr=bvert;

  PetscReal d;
  PointFromPlane(rr,rho1,rho2,rho3,&d);
  if (PetscAbsReal(d)<D_EPS) return(0);

  /* calculate rho_j */
  my_daxpy(ND,-1.0,rr,1,rho1,1);
  my_daxpy(ND,-1.0,rr,1,rho2,1);
  my_daxpy(ND,-1.0,rr,1,rho3,1);

  /* zetal gives ("normal") distance of R from the plane of the triangle */
  zetal=my_ddot(ND,zeta,1,rho1,1);

  /* skip the rest if zetal==0 (R in plane of the triangle)
     -> omega==0 and zetal==0 -> matrix entry=0
  */
  if (PetscAbsReal(zetal)<=D_EPS) {
    return(0);
  }

  rho1l=my_dnrm2(ND,rho1,1);
  rho2l=my_dnrm2(ND,rho2,1);
  rho3l=my_dnrm2(ND,rho3,1);

  PetscReal t_nom,t_denom;
  t_nom=
    rho1l*rho2l*rho3l+
    rho1l*my_ddot(ND,rho2,1,rho3,1)+
    rho2l*my_ddot(ND,rho3,1,rho1,1)+
    rho3l*my_ddot(ND,rho1,1,rho2,1);
  t_denom=
    sqrt(2.0*
      (rho2l*rho3l+my_ddot(ND,rho2,1,rho3,1))*
      (rho3l*rho1l+my_ddot(ND,rho3,1,rho1,1))*
      (rho1l*rho2l+my_ddot(ND,rho1,1,rho2,1))
    );

  /* catch special cases where the argument of acos
     is close to -1.0 or 1.0 or even outside this interval

     use 0.0 instead of D_EPS?
     fixes problems with demag field calculation
     suggested by Hiroki Kobayashi, Fujitsu
  */
  if (t_nom/t_denom<-1.0) {
    omega=(zetal >= 0.0 ? 1.0 : -1.0)*2.0*M_PI;
  }
  /* this case should not occur, but does - e.g. examples1/headfield */
  else if (t_nom/t_denom>1.0) {
    return(0);
  }
  else {
    omega=(zetal >= 0.0 ? 1.0 : -1.0)*2.0*acos(t_nom/t_denom);
  }

  eta1l=my_ddot(ND,eta1,1,rho1,1);
  eta2l=my_ddot(ND,eta2,1,rho2,1);
  eta3l=my_ddot(ND,eta3,1,rho3,1);

  p[0]=log((rho1l+rho2l+s1)/(rho1l+rho2l-s1));
  p[1]=log((rho2l+rho3l+s2)/(rho2l+rho3l-s2));
  p[2]=log((rho3l+rho1l+s3)/(rho3l+rho1l-s3));

  matele[0]=(eta2l*omega-zetal*my_ddot(ND,gamma1,1,p,1))*s2/(8.0*PETSC_PI*a);
  matele[1]=(eta3l*omega-zetal*my_ddot(ND,gamma2,1,p,1))*s3/(8.0*PETSC_PI*a);
  matele[2]=(eta1l*omega-zetal*my_ddot(ND,gamma3,1,p,1))*s1/(8.0*PETSC_PI*a);

  return(0);
}