Please fill out the following information to help in answering your question, and also see tips for posting code snippets. If you don’t provide this information it will take more time to help with your problem!
Geant4 Version: 10.07
Operating System: Window 11
Compiler/Version:
CMake Version:
Dear all,
I am trying to simulate an electric field in Geant4 that uses COMSOL output.
-
I calculated the electric field map using COMSOL with geometry like the attached images and exported them to text files.
-
I constructed my application based on the filed02 and purging_magnet examples. I wrote the application successfully. When I plotted the 3D image using QT, I saw that electric vectors concentrate at a part of the geometry as shown in the attached image. I used the following commands to plot them.
# Draw magneticfield
/vis/scene/add/magneticField
/vis/scene/activateModel Field false
#/vis/scene/add/electricField 20 fullArrow lightArrow
/vis/scene/add/electricField 50 lightArrow
Can anyone help me to solve it? Thank you in advance.
#include "ElectricFieldSetup.hh"
#include "G4SystemOfUnits.hh"
#include "G4AutoLock.hh"
namespace
{
G4Mutex electricField3DLock = G4MUTEX_INITIALIZER;
}
ElectricFieldSetup::ElectricFieldSetup(const char* filename,
double xOffset, double yOffset, double zOffset)
: fXoffset(xOffset), fYoffset(yOffset), fZoffset(zOffset),
invertX(false), invertY(false), invertZ(false)
{
double lenUnit = m;
double fieldUnit = volt/m;
G4cout << "\n-----------------------------------------------------------"
<< "\n Magnetic field"
<< "\n-----------------------------------------------------------";
G4cout << "\n ---> " "Reading the field grid from " << filename << " ... " << G4endl;
/// This is a thread-local class and we have to avoid that all workers open the file at the same time
G4AutoLock lock(&electricField3DLock);
ifstream file(filename); // Open the file for reading.
if (!file.is_open())
{
G4ExceptionDescription ed;
ed << "Could not open input file " << filename << std::endl;
G4Exception("ElectricFieldSetup::ElectricFieldSetup", "electric001", FatalException, ed);
}
// Ignore first blank line
char buffer[256];
file.getline(buffer, 256);
// Read table dimensions
file >> nx >> ny >> nz; // Note dodgy order
G4cout << " [ Number of values x,y,z: "
<< nx << " " << ny << " " << nz << " ] "
<< G4endl;
// Set up storage space for table
xField.resize(nx);
yField.resize(nx);
zField.resize(nx);
int ix, iy, iz;
for(ix=0; ix<nx; ix++)
{
xField[ix].resize(ny);
yField[ix].resize(ny);
zField[ix].resize(ny);
for (iy=0; iy<ny; iy++)
{
xField[ix][iy].resize(nz);
yField[ix][iy].resize(nz);
zField[ix][iy].resize(nz);
}
}
// Ignore other header information
// The first line whose second character is '0' is considered to be the last line of the header.
do
{
file.getline(buffer,256);
}while(buffer[1]!='0');
// Read in the data
double xval, yval, zval, bx, by, bz;
for(ix=0; ix<nx; ix++)
{
for(iy=0; iy<ny; iy++)
{
for (iz=0; iz<nz; iz++)
{
file >> xval >> yval >> zval >> bx >> by >> bz;
//G4cout << " xval: " << xval << " yval: " << yval << " zval: " << zval << endl;
if( ix==0 && iy==0 && iz==0 )
{
minx = xval * lenUnit;
miny = yval * lenUnit;
minz = zval * lenUnit;
}
xField[ix][iy][iz] = bx * fieldUnit;
yField[ix][iy][iz] = by * fieldUnit;
zField[ix][iy][iz] = bz * fieldUnit;
// G4cout << " ix: " << ix
// << " iy: " << iy
// << " iz: " << iz
// << " xField: " << xField[ix][iy][iz]
// << " yField: " << yField[ix][iy][iz]
// << " zField: " << zField[ix][iy][iz]
// << G4endl;
}
}
}
file.close();
lock.unlock();
maxx = xval * lenUnit;
maxy = yval * lenUnit;
maxz = zval * lenUnit;
G4cout << "\n ---> ... done reading " << G4endl;
G4cout << " ---> assumed the order: x, y, z, Bx, By, Bz "
<< "\n ---> Min values x,y,z: " << minx/cm << " " << miny/cm << " " << minz/cm << " cm "
<< "\n ---> Max values x,y,z: " << maxx/cm << " " << maxy/cm << " " << maxz/cm << " cm "
<< "\n ---> The field will be offset by " << zOffset/cm << " cm "
<< G4endl;
// Should really check that the limits are not the wrong way around.
if(maxx < minx) {swap(maxx,minx); invertX = true;}
if(maxy < miny) {swap(maxy,miny); invertY = true;}
if(maxz < minz) {swap(maxz,minz); invertZ = true;}
G4cout << "\nAfter reordering if neccesary"
<< "\n ---> Min values x,y,z: " << minx/mm << " " << miny/mm << " " << minz/mm << " mm "
<< "\n ---> Max values x,y,z: " << maxx/mm << " " << maxy/mm << " " << maxz/mm << " mm ";
dx = maxx - minx;
dy = maxy - miny;
dz = maxz - minz;
G4cout << "\n ---> Dif values x,y,z (range): "
<< dx/cm << " " << dy/cm << " " << dz/cm << " cm in z "
<< "\n-----------------------------------------------------------" << G4endl;
}
void ElectricFieldSetup::GetFieldValue(const double Point[4], double* Bfield) const
{
double x = Point[0] + fXoffset;;
double y = Point[1] + fYoffset;
double z = Point[2] + fZoffset;
// Check that the point is within the defined region
if (x>=minx && x<=maxx &&
y>=miny && y<=maxy &&
z>=minz && z<=maxz)
{
// Position of given point within region, normalized to the range [0,1]
double xfraction = (x - minx) / dx;
double yfraction = (y - miny) / dy;
double zfraction = (z - minz) / dz;
if (invertX) { xfraction = 1 - xfraction;}
if (invertY) { yfraction = 1 - yfraction;}
if (invertZ) { zfraction = 1 - zfraction;}
// Need addresses of these to pass to modf below. modf uses its second argument as an OUTPUT argument.
double xdindex, ydindex, zdindex;
// Position of the point within the cuboid defined by the nearest surrounding tabulated points
double xlocal = ( std::modf(xfraction*(nx-1), &xdindex));
double ylocal = ( std::modf(yfraction*(ny-1), &ydindex));
double zlocal = ( std::modf(zfraction*(nz-1), &zdindex));
// The indices of the nearest tabulated point whose coordinates are all less than those of the given point
int xindex = static_cast<int>(xdindex);
int yindex = static_cast<int>(ydindex);
int zindex = static_cast<int>(zdindex);
#ifdef DEBUG_INTERPOLATING_FIELD
G4cout << "Local x,y,z: " << xlocal << " " << ylocal << " " << zlocal << G4endl;
G4cout << "Index x,y,z: " << xindex << " " << yindex << " " << zindex << G4endl;
double valx0z0, mulx0z0, valx1z0, mulx1z0;
double valx0z1, mulx0z1, valx1z1, mulx1z1;
valx0z0= table[xindex ][0][zindex]; mulx0z0= (1-xlocal) * (1-zlocal);
valx1z0= table[xindex+1][0][zindex]; mulx1z0= xlocal * (1-zlocal);
valx0z1= table[xindex ][0][zindex+1]; mulx0z1= (1-xlocal) * zlocal;
valx1z1= table[xindex+1][0][zindex+1]; mulx1z1= xlocal * zlocal;
#endif
// Full 3-dimensional version
Bfield[0] = 0.0;
Bfield[1] = 0.0;
Bfield[2] = 0.0;
Bfield[3] =
xField[xindex ][yindex ][zindex ] * (1-xlocal) * (1-ylocal) * (1-zlocal) +
xField[xindex ][yindex ][zindex+1] * (1-xlocal) * (1-ylocal) * zlocal +
xField[xindex ][yindex+1][zindex ] * (1-xlocal) * ylocal * (1-zlocal) +
xField[xindex ][yindex+1][zindex+1] * (1-xlocal) * ylocal * zlocal +
xField[xindex+1][yindex ][zindex ] * xlocal * (1-ylocal) * (1-zlocal) +
xField[xindex+1][yindex ][zindex+1] * xlocal * (1-ylocal) * zlocal +
xField[xindex+1][yindex+1][zindex ] * xlocal * ylocal * (1-zlocal) +
xField[xindex+1][yindex+1][zindex+1] * xlocal * ylocal * zlocal ;
Bfield[4] =
yField[xindex ][yindex ][zindex ] * (1-xlocal) * (1-ylocal) * (1-zlocal) +
yField[xindex ][yindex ][zindex+1] * (1-xlocal) * (1-ylocal) * zlocal +
yField[xindex ][yindex+1][zindex ] * (1-xlocal) * ylocal * (1-zlocal) +
yField[xindex ][yindex+1][zindex+1] * (1-xlocal) * ylocal * zlocal +
yField[xindex+1][yindex ][zindex ] * xlocal * (1-ylocal) * (1-zlocal) +
yField[xindex+1][yindex ][zindex+1] * xlocal * (1-ylocal) * zlocal +
yField[xindex+1][yindex+1][zindex ] * xlocal * ylocal * (1-zlocal) +
yField[xindex+1][yindex+1][zindex+1] * xlocal * ylocal * zlocal ;
Bfield[5] =
zField[xindex ][yindex ][zindex ] * (1-xlocal) * (1-ylocal) * (1-zlocal) +
zField[xindex ][yindex ][zindex+1] * (1-xlocal) * (1-ylocal) * zlocal +
zField[xindex ][yindex+1][zindex ] * (1-xlocal) * ylocal * (1-zlocal) +
zField[xindex ][yindex+1][zindex+1] * (1-xlocal) * ylocal * zlocal +
zField[xindex+1][yindex ][zindex ] * xlocal * (1-ylocal) * (1-zlocal) +
zField[xindex+1][yindex ][zindex+1] * xlocal * (1-ylocal) * zlocal +
zField[xindex+1][yindex+1][zindex ] * xlocal * ylocal * (1-zlocal) +
zField[xindex+1][yindex+1][zindex+1] * xlocal * ylocal * zlocal ;
}
else
{
Bfield[3] = 0.0;
Bfield[4] = 0.0;
Bfield[5] = 0.0;
}
}
