yy_mel_util.cc

/* FILE: yy_mel_util.cc                 -*-Mode: c++-*-
 *
 * OOMMF magnetoelastic coupling extension module.
 * yy_mel_util.* contain YY_MELField class , which is to be used by other
 * YY_*MEL classes
 *
 * Release ver. 1.0.1 (2015-03-03)
 * 
 */

#include "yy_mel_util.h"

#include "mesh.h"
#include "rectangularmesh.h"
#include "meshvalue.h"
#include "simstate.h"
#include "threevector.h"
#include "vectorfield.h"
#include "util.h"

#ifdef YY_DEBUG
#include <iostream>
#endif

/* End includes */

YY_MELField::YY_MELField():
  displacement_valid(0), strain_valid(0), 
  MELCoef1_valid(0), MELCoef2_valid(0), pE_pt_buf(0)
{
}

void YY_MELField::Release()
{
  u_cache.Release();
  e_diag_cache.Release();
  e_offdiag_cache.Release();
  u.Release();
  e_diag.Release();
  e_offdiag.Release();
  du.Release();
  de_diag.Release();
  de_offdiag.Release();
  MELCoef1.Release();
  MELCoef2.Release();
}

void YY_MELField::SetMELCoef(const Oxs_SimState& state,
    const Oxs_OwnedPointer<Oxs_ScalarField>& MELCoef1_init,
    const Oxs_OwnedPointer<Oxs_ScalarField>& MELCoef2_init)
{
  MELCoef1_init->FillMeshValue(state.mesh,MELCoef1);
  MELCoef2_init->FillMeshValue(state.mesh,MELCoef2);
  MELCoef1_valid = 1;
  MELCoef2_valid = 1;
}

void YY_MELField::SetDisplacement(const Oxs_SimState& state,
    const Oxs_OwnedPointer<Oxs_VectorField>& u_init)
{
  // Check mesh size
  const OC_INDEX size = state.mesh->Size();
  if(size<1) return;

  u_init->FillMeshValue(state.mesh,u_cache);
  u = u_cache;

  // Initialize du/dt with zeros.
  const ThreeVector zeros(0,0,0);
  du.AdjustSize(state.mesh);
  du=zeros;
  displacement_valid = 1;

  CalculateStrain(state);
  // Update the strain cache.
  e_diag_cache = e_diag;
  e_offdiag_cache = e_offdiag;
  strain_valid = 1;

#ifdef YY_DEBUG
  DisplayValues(state,4,6,0,2,0,2);
#endif
}

void YY_MELField::SetStrain(const Oxs_SimState& state,
    const Oxs_OwnedPointer<Oxs_VectorField>& e_diag_init,
    const Oxs_OwnedPointer<Oxs_VectorField>& e_offdiag_init)
{
  if(state.mesh->Size()<1) return;

  e_diag_init->FillMeshValue(state.mesh,e_diag_cache);
  e_offdiag_init->FillMeshValue(state.mesh,e_offdiag_cache);
  e_diag = e_diag_cache;
  e_offdiag=e_offdiag_cache;

  // Initialize de/dt with zeros.
  const ThreeVector zeros(0,0,0);
  de_diag.AdjustSize(state.mesh);
  de_offdiag.AdjustSize(state.mesh);
  de_diag=zeros;
  de_offdiag=zeros;

  displacement_valid = 0;
  strain_valid = 1;

#ifdef YY_DEBUG
  DisplayValues(state,4,6,0,2,0,2);
#endif
}

void YY_MELField::CalculateStrain(const Oxs_SimState& state)
{
  const Oxs_MeshValue<OC_REAL8m>& Ms = *(state.Ms);
  const Oxs_MeshValue<ThreeVector>& spin = state.spin;
  const Oxs_RectangularMesh* mesh =
    dynamic_cast<const Oxs_RectangularMesh*>(state.mesh);
  const OC_INDEX xdim = mesh->DimX();
  const OC_INDEX ydim = mesh->DimY();
  const OC_INDEX zdim = mesh->DimZ();
  const OC_INDEX xydim = xdim*ydim;
  const OC_REAL8m idelx = 1./mesh->EdgeLengthX();
  const OC_REAL8m idely = 1./mesh->EdgeLengthX();
  const OC_REAL8m idelz = 1./mesh->EdgeLengthX();

  e_diag.AdjustSize(mesh);
  e_offdiag.AdjustSize(mesh);
  de_diag.AdjustSize(mesh);
  de_offdiag.AdjustSize(mesh);

  // Compute du/dx
  // The following calculation assumes that the displacement is only
  // defined in the region where Ms != 0.
  for(OC_INDEX z=0; z<zdim; z++) {
    for(OC_INDEX y=0; y<ydim; y++) {
      for(OC_INDEX x=0; x<xdim; x++) {
        OC_INDEX i = mesh->Index(x,y,z);  // Get base index
        ThreeVector du_dx, ddu_dx;  // ddu_dx = d(du/dx)/dt
        ThreeVector du_dy, ddu_dy;
        ThreeVector du_dz, ddu_dz;
        if(Ms[i]==0.0) {
          e_diag[i].Set(0.0, 0.0, 0.0);
          e_offdiag[i].Set(0.0, 0.0, 0.0);
          de_diag[i].Set(0.0, 0.0, 0.0);
          de_offdiag[i].Set(0.0, 0.0, 0.0);
          continue;
        }

        if(x<xdim-1 && Ms[i+1]!=0.0) { du_dx  = u[i+1]; ddu_dx  = du[i+1]; }
        else                         { du_dx  = u[i];   ddu_dx  = du[i];   }
        if(x>0 && Ms[i-1]!=0.0)      { du_dx -= u[i-1]; ddu_dx -= du[i-1]; }
        else                         { du_dx -= u[i];   ddu_dx -= du[i];   }
        if(x<xdim-1 && Ms[i+1]!=0.0 && x>0 && Ms[i-1]!=0.0) {
          du_dx *= 0.5*idelx; ddu_dx *= 0.5*idelx;
        } else {
          du_dx *= idelx; ddu_dx *= idelx;
        }

        if(y<ydim-1 && Ms[i+xdim]!=0.0) { du_dy  = u[i+xdim]; ddu_dy  = du[i+xdim]; }
        else                            { du_dy  = u[i];      ddu_dy  = du[i];      }
        if(y>0 && Ms[i-xdim]!=0.0)      { du_dy -= u[i-xdim]; ddu_dy -= du[i-xdim]; }
        else                            { du_dy -= u[i];      ddu_dy -= du[i];      }
        if(y<ydim-1 && Ms[i+xdim]!=0.0 && y>0 && Ms[i-xdim]!=0.0) {
          du_dy *= 0.5*idely; ddu_dy *= 0.5*idely;
        } else {
          du_dy *= idely; ddu_dy *= idely;
        }

        if(z<zdim-1 && Ms[i+xydim]!=0.0) { du_dz  = u[i+xydim]; ddu_dz  = du[i+xydim]; }
        else                             { du_dz  = u[i];       ddu_dz  = du[i];       }
        if(z>0 && Ms[i-xydim]!=0.0)      { du_dz -= u[i-xydim]; ddu_dz -= du[i-xydim]; }
        else                             { du_dz -= u[i];       ddu_dz -= du[i];       }
        if(z<zdim-1 && Ms[i+xydim]!=0.0 && z>0 && Ms[i-xydim]!=0.0) {
          du_dz *= 0.5*idelz; ddu_dz *= 0.5*idelz;
        } else {
          du_dz *= idelz; ddu_dz *= idelz;
        }

        e_diag[i].Set(du_dx.x,du_dy.y,du_dz.z);
        e_offdiag[i].Set(
          0.5*(du_dz.y+du_dy.z),
          0.5*(du_dx.z+du_dz.x),
          0.5*(du_dy.x+du_dx.y)
        );
        de_diag[i].Set(ddu_dx.x,ddu_dy.y,ddu_dz.z);
        de_offdiag[i].Set(
          0.5*(ddu_dz.y+ddu_dy.z),
          0.5*(ddu_dx.z+ddu_dz.x),
          0.5*(ddu_dy.x+ddu_dx.y)
        );
      }
    }
  }
}

void YY_MELField::TransformDisplacement(
    const Oxs_SimState& state,
    ThreeVector& row1, ThreeVector& row2, ThreeVector& row3,
    ThreeVector& drow1, ThreeVector& drow2, ThreeVector& drow3)
{
  const OC_INDEX size = state.mesh->Size();
  if(size<1) return;
  if(!displacement_valid) return;

  for(OC_INDEX i=0; i<size; i++) {
    const ThreeVector& v = u_cache[i];
    u[i].Set(row1*v,row2*v,row3*v);
    du[i].Set(drow1*v,drow2*v,drow3*v);
  }
  CalculateStrain(state);

#ifdef YY_DEBUG
  //std::cerr<<"displacement u_cache, u, du/dt"<<endl;
  //std::cerr<<u_cache[20].x<<"  "<<u[20].x<<"  "<<du[20].x<<endl;
  //std::cerr<<u_cache[20].y<<"  "<<u[20].y<<"  "<<du[20].y<<endl;
  //std::cerr<<u_cache[20].z<<"  "<<u[20].z<<"  "<<du[20].z<<endl;
  //std::cerr<<"strain, de/dt"<<endl;
  //std::cerr<<e_diag[20].x<<" "<<e_offdiag[20].z<<" "<<e_offdiag[20].y<<"  ";
  //std::cerr<<de_diag[20].x<<" "<<de_offdiag[20].z<<" "<<de_offdiag[20].y<<endl;
  //std::cerr<<e_offdiag[20].z<<" "<<e_diag[20].y<<" "<<e_offdiag[20].x<<"  ";
  //std::cerr<<de_offdiag[20].z<<" "<<de_diag[20].y<<" "<<de_offdiag[20].x<<endl;
  //std::cerr<<e_offdiag[20].y<<" "<<e_offdiag[20].x<<" "<<e_diag[20].z<<"  ";
  //std::cerr<<de_offdiag[20].y<<" "<<de_offdiag[20].x<<" "<<de_diag[20].z<<endl;
#endif
}

void YY_MELField::TransformStrain(
    const Oxs_SimState& state,
    ThreeVector& row1, ThreeVector& row2, ThreeVector& row3,
    ThreeVector& drow1, ThreeVector& drow2, ThreeVector& drow3)
{
  const OC_INDEX size = state.mesh->Size();
  if(size<1) return;
  if(!strain_valid) return;

  // Transformation by matrix multiplication T*e*transpose(T), where
  //     / row1 \
  // T = | row2 |
  //     \ row3 /
  // and e is the strain matrix.
  ThreeVector trow1, trow2, trow3;  // Temporary matrix (T*e), row 1 through 3.
  ThreeVector dtrow1, dtrow2, dtrow3;  // Time derivative of the above.
  for(OC_INDEX i=0; i<size; i++) {
    const ThreeVector& vd = e_diag_cache[i];
    const ThreeVector& vo = e_offdiag_cache[i];
    // Calculate T*e
    trow1.Set(
        row1.x*vd.x + row1.y*vo.z + row1.z*vo.y,
        row1.x*vo.z + row1.y*vd.y + row1.z*vo.x,
        row1.x*vo.y + row1.y*vo.x + row1.z*vd.z);
    trow2.Set(
        row2.x*vd.x + row2.y*vo.z + row2.z*vo.y,
        row2.x*vo.z + row2.y*vd.y + row2.z*vo.x,
        row2.x*vo.y + row2.y*vo.x + row2.z*vd.z);
    trow3.Set(
        row3.x*vd.x + row3.y*vo.z + row3.z*vo.y,
        row3.x*vo.z + row3.y*vd.y + row3.z*vo.x,
        row3.x*vo.y + row3.y*vo.x + row3.z*vd.z);

    // Calculate (T*e)*transpose(T)
    e_diag[i].Set(
        trow1*row1,
        trow2*row2,
        trow3*row3);
    e_offdiag[i].Set(
        trow2*row3,
        trow1*row3,
        trow1*row2);

    // Calculate d(T*e)/dt
    dtrow1.Set(
        drow1.x*vd.x + drow1.y*vo.z + drow1.z*vo.y,
        drow1.x*vo.z + drow1.y*vd.y + drow1.z*vo.x,
        drow1.x*vo.y + drow1.y*vo.x + drow1.z*vd.z);
    dtrow2.Set(
        drow2.x*vd.x + drow2.y*vo.z + drow2.z*vo.y,
        drow2.x*vo.z + drow2.y*vd.y + drow2.z*vo.x,
        drow2.x*vo.y + drow2.y*vo.x + drow2.z*vd.z);
    dtrow3.Set(
        drow3.x*vd.x + drow3.y*vo.z + drow3.z*vo.y,
        drow3.x*vo.z + drow3.y*vd.y + drow3.z*vo.x,
        drow3.x*vo.y + drow3.y*vo.x + drow3.z*vd.z);

    // Calculate d((T*e)*transpose(T))/dt
    de_diag[i].Set(
        dtrow1*row1 + trow1*drow1,
        dtrow2*row2 + trow2*drow2,
        dtrow3*row3 + trow3*drow3);
    de_offdiag[i].Set(
        dtrow2*row3 + trow2*drow3,
        dtrow1*row3 + trow1*drow3,
        dtrow1*row2 + trow1*drow2);
  }

#ifdef YY_DEBUG
  //std::cerr<<"strain, de/dt"<<endl;
  //std::cerr<<e_diag_cache[20].x<<" "<<e_offdiag_cache[20].z<<" "<<e_offdiag_cache[20].y<<"  ";
  //std::cerr<<e_diag[20].x<<" "<<e_offdiag[20].z<<" "<<e_offdiag[20].y<<"  ";
  //std::cerr<<de_diag[20].x<<" "<<de_offdiag[20].z<<" "<<de_offdiag[20].y<<endl;
  //std::cerr<<e_offdiag_cache[20].z<<" "<<e_diag_cache[20].y<<" "<<e_offdiag_cache[20].x<<"  ";
  //std::cerr<<e_offdiag[20].z<<" "<<e_diag[20].y<<" "<<e_offdiag[20].x<<"  ";
  //std::cerr<<de_offdiag[20].z<<" "<<de_diag[20].y<<" "<<de_offdiag[20].x<<endl;
  //std::cerr<<e_offdiag_cache[20].y<<" "<<e_offdiag_cache[20].x<<" "<<e_diag_cache[20].z<<"  ";
  //std::cerr<<e_offdiag[20].y<<" "<<e_offdiag[20].x<<" "<<e_diag[20].z<<"  ";
  //std::cerr<<de_offdiag[20].y<<" "<<de_offdiag[20].x<<" "<<de_diag[20].z<<endl;
#endif
}

void YY_MELField::Interpolate(
  const Oxs_SimState& state,
  OC_REAL8m working_stage_stopping_time,
  YY_MELField& MELField1,
  YY_MELField& MELField2)
{
  OC_REAL8m coef0, coef1, coef2, coef3;
  OC_REAL8m dcoef0, dcoef1, dcoef2, dcoef3;
  ThreeVector row11, row12, row13;
  ThreeVector drow11, drow12, drow13;
  ThreeVector row21, row22, row23;
  ThreeVector drow21, drow22, drow23;

  OC_REAL8m t = state.stage_elapsed_time / working_stage_stopping_time;
  coef1 = (1-t);
  coef2 = t;
  dcoef1 = -1/working_stage_stopping_time;
  dcoef2 = 1/working_stage_stopping_time;
  row11.Set(coef1,0,0);
  row12.Set(0,coef1,0);
  row13.Set(0,0,coef1);
  drow11.Set(dcoef1,0,0);
  drow12.Set(0,dcoef1,0);
  drow13.Set(0,0,dcoef1);
  row21.Set(coef2,0,0);
  row22.Set(0,coef2,0);
  row23.Set(0,0,coef2);
  drow21.Set(dcoef2,0,0);
  drow22.Set(0,dcoef2,0);
  drow23.Set(0,0,dcoef2);

  MELCoef1 = MELField1.MELCoef1;
  MELCoef2 = MELField1.MELCoef2;
  MELCoef1_valid = MELField1.MELCoef1_valid;
  MELCoef2_valid = MELField1.MELCoef2_valid;
  max_field = MELField1.max_field;
  pE_pt_buf = MELField1.pE_pt_buf;
  if(MELField1.displacement_valid) {
    MELField1.TransformDisplacement(state, row11, row12, row13, drow11, drow12, drow13);
    MELField2.TransformDisplacement(state, row21, row22, row23, drow21, drow22, drow23);
    strain_valid = false;
    u        = MELField1.u;
    u       += MELField2.u;
    u_cache  = MELField1.u_cache;
    u_cache += MELField2.u_cache;
    du       = MELField1.du;
    du      += MELField2.du;
    CalculateStrain(state);
    strain_valid = true;
  } else if(MELField1.strain_valid) {
    MELField1.TransformStrain(state, row11, row12, row13, drow11, drow12, drow13);
    MELField2.TransformStrain(state, row21, row22, row23, drow21, drow22, drow23);
    displacement_valid = false;
    e_diag           = MELField1.e_diag;
    e_diag          += MELField2.e_diag;
    e_diag_cache     = MELField1.e_diag_cache;
    e_diag_cache    += MELField2.e_diag_cache;
    de_diag          = MELField1.de_diag;
    de_diag         += MELField2.de_diag;
    e_offdiag        = MELField1.e_offdiag;
    e_offdiag       += MELField2.e_offdiag;
    e_offdiag_cache = MELField1.e_offdiag_cache;
    e_offdiag_cache += MELField2.e_offdiag_cache;
    de_offdiag       = MELField1.de_offdiag;
    de_offdiag      += MELField2.de_offdiag;
  }
}

void YY_MELField::CalculateMELField(
  const Oxs_SimState& state,
  OC_REAL8m hmult,
  Oxs_MeshValue<ThreeVector>& field_buf,
  Oxs_MeshValue<OC_REAL8m>& energy_buf) const
{
  const OC_INDEX size = state.mesh->Size();
  if(size<1) return;

  const Oxs_MeshValue<ThreeVector>& spin = state.spin;
  const Oxs_MeshValue<OC_REAL8m>& Ms = *(state.Ms);
  const Oxs_MeshValue<OC_REAL8m>& Msi = *(state.Ms_inverse);
  const Oxs_RectangularMesh* mesh =
    dynamic_cast<const Oxs_RectangularMesh*>(state.mesh);
  const OC_INDEX xdim = mesh->DimX();
  const OC_INDEX ydim = mesh->DimX();
  const OC_INDEX zdim = mesh->DimX();
  const OC_INDEX xydim = xdim*ydim;

  // Compute MEL field and energy
  Oxs_MeshValue<ThreeVector> temp_field;
  temp_field.AdjustSize(mesh);
  for(OC_INDEX i=0; i<size; i++) {
    // field_buf[i]*e_diag[i] returns a dot product. Don't use it.
    field_buf[i].x  = spin[i].x*e_diag[i].x;
    field_buf[i].y  = spin[i].y*e_diag[i].y;
    field_buf[i].z  = spin[i].z*e_diag[i].z;
    field_buf[i] *= -1/MU0*2*Msi[i]*MELCoef1[i];
    temp_field[i].x = spin[i].y*e_offdiag[i].z+spin[i].z*e_offdiag[i].y;
    temp_field[i].y = spin[i].x*e_offdiag[i].z+spin[i].z*e_offdiag[i].x;
    temp_field[i].z = spin[i].x*e_offdiag[i].y+spin[i].y*e_offdiag[i].x;
    temp_field[i] *= -1/MU0*2*Msi[i]*MELCoef2[i];
    field_buf[i] += temp_field[i];
  }
  if(hmult != 1.0) field_buf *= hmult;

  // H-field
  max_field.Set(0.,0.,0.);
  if(size>0) {
    OC_INDEX max_i = 0;
    OC_REAL8m max_magsq = field_buf[OC_INDEX(0)].MagSq();
    for(OC_INDEX i=1; i<size; i++) {
      OC_REAL8m magsq = field_buf[i].MagSq();
      if(magsq>max_magsq) {
        max_magsq = magsq;
        max_i = i;
      }
    }
    max_field = field_buf[max_i];
  }

  // Energy density
  OC_REAL8m pE_pt_sum=0;
  for(OC_INDEX i=0; i<size; i++) {
    ThreeVector dH, dH2;
    dH.x  = spin[i].x*de_diag[i].x;
    dH.y  = spin[i].y*de_diag[i].y;
    dH.z  = spin[i].z*de_diag[i].z;
    dH *= -1/MU0*2*Msi[i]*MELCoef1[i];
    dH2.x = spin[i].y*de_offdiag[i].z+spin[i].z*de_offdiag[i].y;
    dH2.y = spin[i].x*de_offdiag[i].z+spin[i].z*de_offdiag[i].x;
    dH2.z = spin[i].x*de_offdiag[i].y+spin[i].y*de_offdiag[i].x;
    dH2 *= -1/MU0*2*Msi[i]*MELCoef2[i];
    dH += dH2;
    ThreeVector temp = (-0.5*MU0*Ms[i])*spin[i];
    energy_buf[i] = field_buf[i]*temp;
    pE_pt_sum += mesh->Volume(i)*(dH*temp);
  }
#ifdef YY_DEBUG
  std::cerr<<"pE_pt: "<<pE_pt_sum<<endl;
#endif
  pE_pt_buf = pE_pt_sum;
}

#ifdef YY_DEBUG
void YY_MELField::DisplayValues(
    const Oxs_SimState& state,
    OC_INDEX xmin, OC_INDEX xmax,
    OC_INDEX ymin, OC_INDEX ymax,
    OC_INDEX zmin, OC_INDEX zmax) const
{
  const Oxs_RectangularMesh* mesh =
    static_cast<const Oxs_RectangularMesh*>(state.mesh);
  const Oxs_MeshValue<OC_REAL8m>& Ms = *(state.Ms);

  // Indices
  fprintf(stderr,"Indices:\n");
  for(OC_INDEX y=ymin; y<ymax+1; y++) {
    for(OC_INDEX z=zmin; z<zmax+1; z++) {
      for(OC_INDEX x=xmin; x<xmax+1; x++) {
        OC_INDEX i = mesh->Index(x,y,z);
        fprintf(stderr,"%ld %ld %ld %ld ",x,y,z,i);
      }
      fprintf(stderr,"| ");
    }
    fprintf(stderr,"\n");
  }
  fprintf(stderr,"\n");

  // Ms
  fprintf(stderr,"Ms:\n");
  for(OC_INDEX y=ymin; y<ymax+1; y++) {
    for(OC_INDEX z=zmin; z<zmax+1; z++) {
      for(OC_INDEX x=xmin; x<xmax+1; x++) {
        OC_INDEX i = mesh->Index(x,y,z);
        fprintf(stderr,"%e ",Ms[i]);
      }
      fprintf(stderr,"| ");
    }
    fprintf(stderr,"\n");
  }
  fprintf(stderr,"\n");

  if(displacement_valid) {
    // u.x
    fprintf(stderr,"u.x:\n");
    for(OC_INDEX y=ymin; y<ymax+1; y++) {
      for(OC_INDEX z=zmin; z<zmax+1; z++) {
        for(OC_INDEX x=xmin; x<xmax+1; x++) {
          OC_INDEX i = mesh->Index(x,y,z);
          fprintf(stderr,"%e ",u[i].x);
        }
        fprintf(stderr,"| ");
      }
      fprintf(stderr,"\n");
    }
    fprintf(stderr,"\n");

    // u.y
    fprintf(stderr,"u.y:\n");
    for(OC_INDEX y=ymin; y<ymax+1; y++) {
      for(OC_INDEX z=zmin; z<zmax+1; z++) {
        for(OC_INDEX x=xmin; x<xmax+1; x++) {
          OC_INDEX i = mesh->Index(x,y,z);
          fprintf(stderr,"%e ",u[i].y);
        }
        fprintf(stderr,"| ");
      }
      fprintf(stderr,"\n");
    }
    fprintf(stderr,"\n");

    // u.z
    fprintf(stderr,"u.z:\n");
    for(OC_INDEX y=ymin; y<ymax+1; y++) {
      for(OC_INDEX z=zmin; z<zmax+1; z++) {
        for(OC_INDEX x=xmin; x<xmax+1; x++) {
          OC_INDEX i = mesh->Index(x,y,z);
          fprintf(stderr,"%e ",u[i].z);
        }
        fprintf(stderr,"| ");
      }
      fprintf(stderr,"\n");
    }
    fprintf(stderr,"\n");
  }

  // e_diag.x
  fprintf(stderr,"e_diag.x:\n");
  for(OC_INDEX y=ymin; y<ymax+1; y++) {
    for(OC_INDEX z=zmin; z<zmax+1; z++) {
      for(OC_INDEX x=xmin; x<xmax+1; x++) {
        OC_INDEX i = mesh->Index(x,y,z);
        fprintf(stderr,"%e ",e_diag[i].x);
      }
      fprintf(stderr,"| ");
    }
    fprintf(stderr,"\n");
  }
  fprintf(stderr,"\n");

  // e_diag.y
  fprintf(stderr,"e_diag.y:\n");
  for(OC_INDEX y=ymin; y<ymax+1; y++) {
    for(OC_INDEX z=zmin; z<zmax+1; z++) {
      for(OC_INDEX x=xmin; x<xmax+1; x++) {
        OC_INDEX i = mesh->Index(x,y,z);
        fprintf(stderr,"%e ",e_diag[i].y);
      }
      fprintf(stderr,"| ");
    }
    fprintf(stderr,"\n");
  }
  fprintf(stderr,"\n");

  // e_diag.z
  fprintf(stderr,"e_diag.z:\n");
  for(OC_INDEX y=ymin; y<ymax+1; y++) {
    for(OC_INDEX z=zmin; z<zmax+1; z++) {
      for(OC_INDEX x=xmin; x<xmax+1; x++) {
        OC_INDEX i = mesh->Index(x,y,z);
        fprintf(stderr,"%e ",e_diag[i].z);
      }
      fprintf(stderr,"| ");
    }
    fprintf(stderr,"\n");
  }
  fprintf(stderr,"\n");

  // e_offdiag.x
  fprintf(stderr,"e_offdiag.x:\n");
  for(OC_INDEX y=ymin; y<ymax+1; y++) {
    for(OC_INDEX z=zmin; z<zmax+1; z++) {
      for(OC_INDEX x=xmin; x<xmax+1; x++) {
        OC_INDEX i = mesh->Index(x,y,z);
        fprintf(stderr,"%e ",e_offdiag[i].x);
      }
      fprintf(stderr,"| ");
    }
    fprintf(stderr,"\n");
  }
  fprintf(stderr,"\n");

  // e_offdiag.y
  fprintf(stderr,"e_offdiag.y:\n");
  for(OC_INDEX y=ymin; y<ymax+1; y++) {
    for(OC_INDEX z=zmin; z<zmax+1; z++) {
      for(OC_INDEX x=xmin; x<xmax+1; x++) {
        OC_INDEX i = mesh->Index(x,y,z);
        fprintf(stderr,"%e ",e_offdiag[i].y);
      }
      fprintf(stderr,"| ");
    }
    fprintf(stderr,"\n");
  }
  fprintf(stderr,"\n");

  // e_offdiag.z
  fprintf(stderr,"e_offdiag.z:\n");
  for(OC_INDEX y=ymin; y<ymax+1; y++) {
    for(OC_INDEX z=zmin; z<zmax+1; z++) {
      for(OC_INDEX x=xmin; x<xmax+1; x++) {
        OC_INDEX i = mesh->Index(x,y,z);
        fprintf(stderr,"%e ",e_offdiag[i].z);
      }
      fprintf(stderr,"| ");
    }
    fprintf(stderr,"\n");
  }
  fprintf(stderr,"\n");

}
#endif  // YY_DEBUG