Garfield++ v1r0
A toolkit for the detailed simulation of particle detectors based on ionisation measurement in gases and semiconductors
Loading...
Searching...
No Matches
Classes | Public Member Functions | Protected Member Functions | Protected Attributes | Friends | List of all members
Garfield::ComponentFieldMap Class Referenceabstract

#include <ComponentFieldMap.hh>

+ Inheritance diagram for Garfield::ComponentFieldMap:

Classes

struct  element
 
struct  material
 
struct  node
 

Public Member Functions

 ComponentFieldMap ()
 
virtual ~ComponentFieldMap ()
 
virtual void SetRange ()
 
void PrintRange ()
 
virtual bool IsInBoundingBox (const double x, const double y, const double z)
 
virtual bool GetBoundingBox (double &xmin, double &ymin, double &zmin, double &xmax, double &ymax, double &zmax)
 
bool GetVoltageRange (double &vmin, double &vmax)
 
void PrintMaterials ()
 
void DriftMedium (int imat)
 
void NotDriftMedium (int imat)
 
int GetNumberOfMaterials ()
 
double GetPermittivity (const int imat)
 
double GetConductivity (const int imat)
 
void SetMedium (const int imat, Medium *medium)
 
Medium * GetMedium (const unsigned int &i) const
 
Medium * GetMedium (const double &x, const double &y, const double &z)=0
 
int GetNumberOfMedia ()
 
int GetNumberOfElements () const
 
bool GetElement (const int i, double &vol, double &dmin, double &dmax)
 
virtual void ElectricField (const double x, const double y, const double z, double &ex, double &ey, double &ez, Medium *&m, int &status)=0
 
virtual void ElectricField (const double x, const double y, const double z, double &ex, double &ey, double &ez, double &v, Medium *&m, int &status)=0
 
virtual void WeightingField (const double x, const double y, const double z, double &wx, double &wy, double &wz, const std::string label)=0
 
virtual double WeightingPotential (const double x, const double y, const double z, const std::string label)=0
 
void EnableCheckMapIndices ()
 
void DisableCheckMapIndices ()
 
void EnableDeleteBackgroundElements ()
 
void DisableDeleteBackgroundElements ()
 
- Public Member Functions inherited from Garfield::ComponentBase
 ComponentBase ()
 
virtual ~ComponentBase ()
 
virtual void SetGeometry (GeometryBase *geo)
 
virtual void Clear ()
 
virtual Medium * GetMedium (const double &x, const double &y, const double &z)
 
virtual void ElectricField (const double x, const double y, const double z, double &ex, double &ey, double &ez, Medium *&m, int &status)=0
 
virtual void ElectricField (const double x, const double y, const double z, double &ex, double &ey, double &ez, double &v, Medium *&m, int &status)=0
 
virtual bool GetVoltageRange (double &vmin, double &vmax)=0
 
virtual void WeightingField (const double x, const double y, const double z, double &wx, double &wy, double &wz, const std::string label)
 
virtual double WeightingPotential (const double x, const double y, const double z, const std::string label)
 
virtual void MagneticField (const double x, const double y, const double z, double &bx, double &by, double &bz, int &status)
 
void SetMagneticField (const double bx, const double by, const double bz)
 
virtual bool IsReady ()
 
virtual bool GetBoundingBox (double &xmin, double &ymin, double &zmin, double &xmax, double &ymax, double &zmax)
 
virtual bool IsWireCrossed (const double x0, const double y0, const double z0, const double x1, const double y1, const double z1, double &xc, double &yc, double &zc)
 
virtual bool IsInTrapRadius (double x0, double y0, double z0, double &xw, double &yw, double &rw)
 
void EnablePeriodicityX ()
 
void DisablePeriodicityX ()
 
void EnablePeriodicityY ()
 
void DisablePeriodicityY ()
 
void EnablePeriodicityZ ()
 
void DisablePeriodicityZ ()
 
void EnableMirrorPeriodicityX ()
 
void DisableMirrorPeriodicityX ()
 
void EnableMirrorPeriodicityY ()
 
void DisableMirrorPeriodicityY ()
 
void EnableMirrorPeriodicityZ ()
 
void DisableMirrorPeriodicityZ ()
 
void EnableAxialPeriodicityX ()
 
void DisableAxialPeriodicityX ()
 
void EnableAxialPeriodicityY ()
 
void DisableAxialPeriodicityY ()
 
void EnableAxialPeriodicityZ ()
 
void DisableAxialPeriodicityZ ()
 
void EnableRotationSymmetryX ()
 
void DisableRotationSymmetryX ()
 
void EnableRotationSymmetryY ()
 
void DisableRotationSymmetryY ()
 
void EnableRotationSymmetryZ ()
 
void DisableRotationSymmetryZ ()
 
void EnableDebugging ()
 
void DisableDebugging ()
 

Protected Member Functions

void Reset ()
 
virtual void UpdatePeriodicity ()=0
 
void UpdatePeriodicity2d ()
 
void UpdatePeriodicityCommon ()
 
int Coordinates3 (double x, double y, double z, double &t1, double &t2, double &t3, double &t4, double jac[4][4], double &det, int imap)
 
int Coordinates4 (double x, double y, double z, double &t1, double &t2, double &t3, double &t4, double jac[4][4], double &det, int imap)
 
int Coordinates5 (double x, double y, double z, double &t1, double &t2, double &t3, double &t4, double jac[4][4], double &det, int imap)
 
int Coordinates12 (double x, double y, double z, double &t1, double &t2, double &t3, double &t4, int imap)
 
int Coordinates13 (double x, double y, double z, double &t1, double &t2, double &t3, double &t4, double jac[4][4], double &det, int imap)
 
int CoordinatesCube (double x, double y, double z, double &t1, double &t2, double &t3, TMatrixD *&jac, std::vector< TMatrixD * > &dN, int imap)
 
void Jacobian3 (int i, double u, double v, double w, double &det, double jac[4][4])
 
void Jacobian5 (int i, double u, double v, double &det, double jac[4][4])
 
void Jacobian13 (int i, double t, double u, double v, double w, double &det, double jac[4][4])
 
void JacobianCube (int i, double t1, double t2, double t3, TMatrixD *&jac, std::vector< TMatrixD * > &dN)
 
int FindElement5 (const double x, const double y, const double z, double &t1, double &t2, double &t3, double &t4, double jac[4][4], double &det)
 
int FindElement13 (const double x, const double y, const double z, double &t1, double &t2, double &t3, double &t4, double jac[4][4], double &det)
 
int FindElementCube (const double x, const double y, const double z, double &t1, double &t2, double &t3, TMatrixD *&jac, std::vector< TMatrixD * > &dN)
 
void MapCoordinates (double &xpos, double &ypos, double &zpos, bool &xmirrored, bool &ymirrored, bool &zmirrored, double &rcoordinate, double &rotation) const
 
void UnmapFields (double &ex, double &ey, double &ez, double &xpos, double &ypos, double &zpos, bool &xmirrored, bool &ymirrored, bool &zmirrored, double &rcoordinate, double &rotation)
 
int ReadInteger (char *token, int def, bool &error)
 
double ReadDouble (char *token, double def, bool &error)
 
virtual double GetElementVolume (const int i)=0
 
virtual void GetAspectRatio (const int i, double &dmin, double &dmax)=0
 
void CalculateElementBoundingBoxes (void)
 
virtual void Reset ()=0
 
virtual void UpdatePeriodicity ()=0
 

Protected Attributes

bool is3d
 
int nElements
 
std::vector< element > elements
 
int lastElement
 
bool cacheElemBoundingBoxes
 
int nNodes
 
std::vector< node > nodes
 
int nMaterials
 
std::vector< material > materials
 
int nWeightingFields
 
std::vector< std::string > wfields
 
std::vector< bool > wfieldsOk
 
bool hasBoundingBox
 
double xMinBoundingBox
 
double yMinBoundingBox
 
double zMinBoundingBox
 
double xMaxBoundingBox
 
double yMaxBoundingBox
 
double zMaxBoundingBox
 
double mapxmin
 
double mapymin
 
double mapzmin
 
double mapxmax
 
double mapymax
 
double mapzmax
 
double mapxamin
 
double mapyamin
 
double mapzamin
 
double mapxamax
 
double mapyamax
 
double mapzamax
 
double mapvmin
 
double mapvmax
 
bool setangx
 
bool setangy
 
bool setangz
 
double mapsx
 
double mapsy
 
double mapsz
 
double cellsx
 
double cellsy
 
double cellsz
 
double mapnxa
 
double mapnya
 
double mapnza
 
bool deleteBackground
 
bool checkMultipleElement
 
bool warning
 
- Protected Attributes inherited from Garfield::ComponentBase
std::string m_className
 
GeometryBase * theGeometry
 
bool ready
 
bool xPeriodic
 
bool yPeriodic
 
bool zPeriodic
 
bool xMirrorPeriodic
 
bool yMirrorPeriodic
 
bool zMirrorPeriodic
 
bool xAxiallyPeriodic
 
bool yAxiallyPeriodic
 
bool zAxiallyPeriodic
 
bool xRotationSymmetry
 
bool yRotationSymmetry
 
bool zRotationSymmetry
 
double bx0
 
double by0
 
double bz0
 
bool debug
 

Friends

class ViewFEMesh
 

Detailed Description

Definition at line 9 of file ComponentFieldMap.hh.

Constructor & Destructor Documentation

◆ ComponentFieldMap()

Garfield::ComponentFieldMap::ComponentFieldMap ( )

Definition at line 15 of file ComponentFieldMap.cc.

16 : is3d(true),
17 nElements(-1),
18 lastElement(-1),
19 cacheElemBoundingBoxes(false),
20 nNodes(-1),
21 nMaterials(-1),
22 nWeightingFields(0),
23 hasBoundingBox(false),
24 deleteBackground(true),
25 checkMultipleElement(false),
26 warning(false) {
27
28 m_className = "ComponentFieldMap";
29
30 materials.clear();
31 elements.clear();
32 nodes.clear();
33 wfields.clear();
34 wfieldsOk.clear();
35}
Garfield::ComponentFieldMap::materials
std::vector< material > materials
Definition: ComponentFieldMap.hh:124
Garfield::ComponentFieldMap::elements
std::vector< element > elements
Definition: ComponentFieldMap.hh:96
Garfield::ComponentFieldMap::wfields
std::vector< std::string > wfields
Definition: ComponentFieldMap.hh:127

◆ ~ComponentFieldMap()

virtual Garfield::ComponentFieldMap::~ComponentFieldMap ( )
inlinevirtual

Definition at line 15 of file ComponentFieldMap.hh.

15{}

Member Function Documentation

◆ CalculateElementBoundingBoxes()

void Garfield::ComponentFieldMap::CalculateElementBoundingBoxes ( void  )
protected

Definition at line 4249 of file ComponentFieldMap.cc.

4249 {
4250
4251 // Do not proceed if not properly initialised.
4252 if (!ready) {
4253 std::cerr << m_className << "::CalculateElementBoundingBoxes:\n";
4254 std::cerr << " Field map not yet initialised.\n";
4255 std::cerr << " Bounding boxes of elements cannot be computed.\n";
4256 return;
4257 }
4258
4259 // Calculate the bounding boxes of all elements
4260 for (int i = 0; i < nElements; ++i) {
4261 element& elem = elements[i];
4262 elem.xmin = std::min(
4263 std::min(nodes[elem.emap[0]].x, nodes[elem.emap[1]].x),
4264 std::min(nodes[elem.emap[2]].x, nodes[elem.emap[3]].x));
4265 elem.xmax = std::max(
4266 std::max(nodes[elem.emap[0]].x, nodes[elem.emap[1]].x),
4267 std::max(nodes[elem.emap[2]].x, nodes[elem.emap[3]].x));
4268 elem.ymin = std::min(
4269 std::min(nodes[elem.emap[0]].y, nodes[elem.emap[1]].y),
4270 std::min(nodes[elem.emap[2]].y, nodes[elem.emap[3]].y));
4271 elem.ymax = std::max(
4272 std::max(nodes[elem.emap[0]].y, nodes[elem.emap[1]].y),
4273 std::max(nodes[elem.emap[2]].y, nodes[elem.emap[3]].y));
4274 elem.zmin = std::min(
4275 std::min(nodes[elem.emap[0]].z, nodes[elem.emap[1]].z),
4276 std::min(nodes[elem.emap[2]].z, nodes[elem.emap[3]].z));
4277 elem.zmax = std::max(
4278 std::max(nodes[elem.emap[0]].z, nodes[elem.emap[1]].z),
4279 std::max(nodes[elem.emap[2]].z, nodes[elem.emap[3]].z));
4280 }
4281}

Referenced by FindElement13(), and FindElement5().

◆ Coordinates12()

int Garfield::ComponentFieldMap::Coordinates12 ( double  x,
double  y,
double  z,
double &  t1,
double &  t2,
double &  t3,
double &  t4,
int  imap 
)
protected

Definition at line 3005 of file ComponentFieldMap.cc.

3007 {
3008
3009 if (debug) {
3010 std::cout << m_className << "::Coordinates12:\n";
3011 std::cout << " Point (" << x << ", " << y << ", " << z << ") for element "
3012 << imap << "\n";
3013 }
3014
3015 // Failure flag
3016 int ifail = 1;
3017
3018 // Compute tetrahedral coordinates.
3019 t1 =
3020 (x - nodes[elements[imap].emap[1]].x) *
3021 ((nodes[elements[imap].emap[2]].y - nodes[elements[imap].emap[1]].y) *
3022 (nodes[elements[imap].emap[3]].z -
3023 nodes[elements[imap].emap[1]].z) -
3024 (nodes[elements[imap].emap[3]].y - nodes[elements[imap].emap[1]].y) *
3025 (nodes[elements[imap].emap[2]].z -
3026 nodes[elements[imap].emap[1]].z)) +
3027 (y - nodes[elements[imap].emap[1]].y) *
3028 ((nodes[elements[imap].emap[2]].z - nodes[elements[imap].emap[1]].z) *
3029 (nodes[elements[imap].emap[3]].x -
3030 nodes[elements[imap].emap[1]].x) -
3031 (nodes[elements[imap].emap[3]].z - nodes[elements[imap].emap[1]].z) *
3032 (nodes[elements[imap].emap[2]].x -
3033 nodes[elements[imap].emap[1]].x)) +
3034 (z - nodes[elements[imap].emap[1]].z) *
3035 ((nodes[elements[imap].emap[2]].x - nodes[elements[imap].emap[1]].x) *
3036 (nodes[elements[imap].emap[3]].y -
3037 nodes[elements[imap].emap[1]].y) -
3038 (nodes[elements[imap].emap[3]].x - nodes[elements[imap].emap[1]].x) *
3039 (nodes[elements[imap].emap[2]].y -
3040 nodes[elements[imap].emap[1]].y));
3041 t2 =
3042 (x - nodes[elements[imap].emap[2]].x) *
3043 ((nodes[elements[imap].emap[0]].y - nodes[elements[imap].emap[2]].y) *
3044 (nodes[elements[imap].emap[3]].z -
3045 nodes[elements[imap].emap[2]].z) -
3046 (nodes[elements[imap].emap[3]].y - nodes[elements[imap].emap[2]].y) *
3047 (nodes[elements[imap].emap[0]].z -
3048 nodes[elements[imap].emap[2]].z)) +
3049 (y - nodes[elements[imap].emap[2]].y) *
3050 ((nodes[elements[imap].emap[0]].z - nodes[elements[imap].emap[2]].z) *
3051 (nodes[elements[imap].emap[3]].x -
3052 nodes[elements[imap].emap[2]].x) -
3053 (nodes[elements[imap].emap[3]].z - nodes[elements[imap].emap[2]].z) *
3054 (nodes[elements[imap].emap[0]].x -
3055 nodes[elements[imap].emap[2]].x)) +
3056 (z - nodes[elements[imap].emap[2]].z) *
3057 ((nodes[elements[imap].emap[0]].x - nodes[elements[imap].emap[2]].x) *
3058 (nodes[elements[imap].emap[3]].y -
3059 nodes[elements[imap].emap[2]].y) -
3060 (nodes[elements[imap].emap[3]].x - nodes[elements[imap].emap[2]].x) *
3061 (nodes[elements[imap].emap[0]].y -
3062 nodes[elements[imap].emap[2]].y));
3063 t3 =
3064 (x - nodes[elements[imap].emap[3]].x) *
3065 ((nodes[elements[imap].emap[0]].y - nodes[elements[imap].emap[3]].y) *
3066 (nodes[elements[imap].emap[1]].z -
3067 nodes[elements[imap].emap[3]].z) -
3068 (nodes[elements[imap].emap[1]].y - nodes[elements[imap].emap[3]].y) *
3069 (nodes[elements[imap].emap[0]].z -
3070 nodes[elements[imap].emap[3]].z)) +
3071 (y - nodes[elements[imap].emap[3]].y) *
3072 ((nodes[elements[imap].emap[0]].z - nodes[elements[imap].emap[3]].z) *
3073 (nodes[elements[imap].emap[1]].x -
3074 nodes[elements[imap].emap[3]].x) -
3075 (nodes[elements[imap].emap[1]].z - nodes[elements[imap].emap[3]].z) *
3076 (nodes[elements[imap].emap[0]].x -
3077 nodes[elements[imap].emap[3]].x)) +
3078 (z - nodes[elements[imap].emap[3]].z) *
3079 ((nodes[elements[imap].emap[0]].x - nodes[elements[imap].emap[3]].x) *
3080 (nodes[elements[imap].emap[1]].y -
3081 nodes[elements[imap].emap[3]].y) -
3082 (nodes[elements[imap].emap[1]].x - nodes[elements[imap].emap[3]].x) *
3083 (nodes[elements[imap].emap[0]].y -
3084 nodes[elements[imap].emap[3]].y));
3085 t4 =
3086 (x - nodes[elements[imap].emap[0]].x) *
3087 ((nodes[elements[imap].emap[2]].y - nodes[elements[imap].emap[0]].y) *
3088 (nodes[elements[imap].emap[1]].z -
3089 nodes[elements[imap].emap[0]].z) -
3090 (nodes[elements[imap].emap[1]].y - nodes[elements[imap].emap[0]].y) *
3091 (nodes[elements[imap].emap[2]].z -
3092 nodes[elements[imap].emap[0]].z)) +
3093 (y - nodes[elements[imap].emap[0]].y) *
3094 ((nodes[elements[imap].emap[2]].z - nodes[elements[imap].emap[0]].z) *
3095 (nodes[elements[imap].emap[1]].x -
3096 nodes[elements[imap].emap[0]].x) -
3097 (nodes[elements[imap].emap[1]].z - nodes[elements[imap].emap[0]].z) *
3098 (nodes[elements[imap].emap[2]].x -
3099 nodes[elements[imap].emap[0]].x)) +
3100 (z - nodes[elements[imap].emap[0]].z) *
3101 ((nodes[elements[imap].emap[2]].x - nodes[elements[imap].emap[0]].x) *
3102 (nodes[elements[imap].emap[1]].y -
3103 nodes[elements[imap].emap[0]].y) -
3104 (nodes[elements[imap].emap[1]].x - nodes[elements[imap].emap[0]].x) *
3105 (nodes[elements[imap].emap[2]].y -
3106 nodes[elements[imap].emap[0]].y));
3107 t1 = t1 /
3108 ((nodes[elements[imap].emap[0]].x - nodes[elements[imap].emap[1]].x) *
3109 ((nodes[elements[imap].emap[2]].y -
3110 nodes[elements[imap].emap[1]].y) *
3111 (nodes[elements[imap].emap[3]].z -
3112 nodes[elements[imap].emap[1]].z) -
3113 (nodes[elements[imap].emap[3]].y -
3114 nodes[elements[imap].emap[1]].y) *
3115 (nodes[elements[imap].emap[2]].z -
3116 nodes[elements[imap].emap[1]].z)) +
3117 (nodes[elements[imap].emap[0]].y - nodes[elements[imap].emap[1]].y) *
3118 ((nodes[elements[imap].emap[2]].z -
3119 nodes[elements[imap].emap[1]].z) *
3120 (nodes[elements[imap].emap[3]].x -
3121 nodes[elements[imap].emap[1]].x) -
3122 (nodes[elements[imap].emap[3]].z -
3123 nodes[elements[imap].emap[1]].z) *
3124 (nodes[elements[imap].emap[2]].x -
3125 nodes[elements[imap].emap[1]].x)) +
3126 (nodes[elements[imap].emap[0]].z - nodes[elements[imap].emap[1]].z) *
3127 ((nodes[elements[imap].emap[2]].x -
3128 nodes[elements[imap].emap[1]].x) *
3129 (nodes[elements[imap].emap[3]].y -
3130 nodes[elements[imap].emap[1]].y) -
3131 (nodes[elements[imap].emap[3]].x -
3132 nodes[elements[imap].emap[1]].x) *
3133 (nodes[elements[imap].emap[2]].y -
3134 nodes[elements[imap].emap[1]].y)));
3135 t2 = t2 /
3136 ((nodes[elements[imap].emap[1]].x - nodes[elements[imap].emap[2]].x) *
3137 ((nodes[elements[imap].emap[0]].y -
3138 nodes[elements[imap].emap[2]].y) *
3139 (nodes[elements[imap].emap[3]].z -
3140 nodes[elements[imap].emap[2]].z) -
3141 (nodes[elements[imap].emap[3]].y -
3142 nodes[elements[imap].emap[2]].y) *
3143 (nodes[elements[imap].emap[0]].z -
3144 nodes[elements[imap].emap[2]].z)) +
3145 (nodes[elements[imap].emap[1]].y - nodes[elements[imap].emap[2]].y) *
3146 ((nodes[elements[imap].emap[0]].z -
3147 nodes[elements[imap].emap[2]].z) *
3148 (nodes[elements[imap].emap[3]].x -
3149 nodes[elements[imap].emap[2]].x) -
3150 (nodes[elements[imap].emap[3]].z -
3151 nodes[elements[imap].emap[2]].z) *
3152 (nodes[elements[imap].emap[0]].x -
3153 nodes[elements[imap].emap[2]].x)) +
3154 (nodes[elements[imap].emap[1]].z - nodes[elements[imap].emap[2]].z) *
3155 ((nodes[elements[imap].emap[0]].x -
3156 nodes[elements[imap].emap[2]].x) *
3157 (nodes[elements[imap].emap[3]].y -
3158 nodes[elements[imap].emap[2]].y) -
3159 (nodes[elements[imap].emap[3]].x -
3160 nodes[elements[imap].emap[2]].x) *
3161 (nodes[elements[imap].emap[0]].y -
3162 nodes[elements[imap].emap[2]].y)));
3163 t3 = t3 /
3164 ((nodes[elements[imap].emap[2]].x - nodes[elements[imap].emap[3]].x) *
3165 ((nodes[elements[imap].emap[0]].y -
3166 nodes[elements[imap].emap[3]].y) *
3167 (nodes[elements[imap].emap[1]].z -
3168 nodes[elements[imap].emap[3]].z) -
3169 (nodes[elements[imap].emap[1]].y -
3170 nodes[elements[imap].emap[3]].y) *
3171 (nodes[elements[imap].emap[0]].z -
3172 nodes[elements[imap].emap[3]].z)) +
3173 (nodes[elements[imap].emap[2]].y - nodes[elements[imap].emap[3]].y) *
3174 ((nodes[elements[imap].emap[0]].z -
3175 nodes[elements[imap].emap[3]].z) *
3176 (nodes[elements[imap].emap[1]].x -
3177 nodes[elements[imap].emap[3]].x) -
3178 (nodes[elements[imap].emap[1]].z -
3179 nodes[elements[imap].emap[3]].z) *
3180 (nodes[elements[imap].emap[0]].x -
3181 nodes[elements[imap].emap[3]].x)) +
3182 (nodes[elements[imap].emap[2]].z - nodes[elements[imap].emap[3]].z) *
3183 ((nodes[elements[imap].emap[0]].x -
3184 nodes[elements[imap].emap[3]].x) *
3185 (nodes[elements[imap].emap[1]].y -
3186 nodes[elements[imap].emap[3]].y) -
3187 (nodes[elements[imap].emap[1]].x -
3188 nodes[elements[imap].emap[3]].x) *
3189 (nodes[elements[imap].emap[0]].y -
3190 nodes[elements[imap].emap[3]].y)));
3191 t4 = t4 /
3192 ((nodes[elements[imap].emap[3]].x - nodes[elements[imap].emap[0]].x) *
3193 ((nodes[elements[imap].emap[2]].y -
3194 nodes[elements[imap].emap[0]].y) *
3195 (nodes[elements[imap].emap[1]].z -
3196 nodes[elements[imap].emap[0]].z) -
3197 (nodes[elements[imap].emap[1]].y -
3198 nodes[elements[imap].emap[0]].y) *
3199 (nodes[elements[imap].emap[2]].z -
3200 nodes[elements[imap].emap[0]].z)) +
3201 (nodes[elements[imap].emap[3]].y - nodes[elements[imap].emap[0]].y) *
3202 ((nodes[elements[imap].emap[2]].z -
3203 nodes[elements[imap].emap[0]].z) *
3204 (nodes[elements[imap].emap[1]].x -
3205 nodes[elements[imap].emap[0]].x) -
3206 (nodes[elements[imap].emap[1]].z -
3207 nodes[elements[imap].emap[0]].z) *
3208 (nodes[elements[imap].emap[2]].x -
3209 nodes[elements[imap].emap[0]].x)) +
3210 (nodes[elements[imap].emap[3]].z - nodes[elements[imap].emap[0]].z) *
3211 ((nodes[elements[imap].emap[2]].x -
3212 nodes[elements[imap].emap[0]].x) *
3213 (nodes[elements[imap].emap[1]].y -
3214 nodes[elements[imap].emap[0]].y) -
3215 (nodes[elements[imap].emap[1]].x -
3216 nodes[elements[imap].emap[0]].x) *
3217 (nodes[elements[imap].emap[2]].y -
3218 nodes[elements[imap].emap[0]].y)));
3219
3220 // Result
3221 if (debug) {
3222 std::cout << m_className << "::Coordinates12:\n";
3223 std::cout << " Tetrahedral coordinates (t, u, v, w) = (" << t1 << ", "
3224 << t2 << ", " << t3 << ", " << t4
3225 << ") sum = " << t1 + t2 + t3 + t4 << ".\n";
3226 }
3227 // Re-compute the (x,y,z) position for this coordinate.
3228 if (debug) {
3229 double xr = nodes[elements[imap].emap[0]].x * t1 +
3230 nodes[elements[imap].emap[1]].x * t2 +
3231 nodes[elements[imap].emap[2]].x * t3 +
3232 nodes[elements[imap].emap[3]].x * t4;
3233 double yr = nodes[elements[imap].emap[0]].y * t1 +
3234 nodes[elements[imap].emap[1]].y * t2 +
3235 nodes[elements[imap].emap[2]].y * t3 +
3236 nodes[elements[imap].emap[3]].y * t4;
3237 double zr = nodes[elements[imap].emap[0]].z * t1 +
3238 nodes[elements[imap].emap[1]].z * t2 +
3239 nodes[elements[imap].emap[2]].z * t3 +
3240 nodes[elements[imap].emap[3]].z * t4;
3241 double sr = t1 + t2 + t3 + t4;
3242 std::cout << m_className << "::Coordinates12:\n";
3243 std::cout << " Position requested: (" << x << ", " << y << ", " << z
3244 << ")\n";
3245 std::cout << " Reconstructed: (" << xr << ", " << yr << ", "
3246 << zr << ")\n";
3247 std::cout << " Difference: (" << x - xr << ", " << y - yr
3248 << ", " << z - zr << ")\n";
3249 std::cout << " Checksum - 1: " << sr - 1 << "\n";
3250 }
3251
3252 // This should always work.
3253 ifail = 0;
3254 return ifail;
3255}

Referenced by Coordinates13().

◆ Coordinates13()

int Garfield::ComponentFieldMap::Coordinates13 ( double  x,
double  y,
double  z,
double &  t1,
double &  t2,
double &  t3,
double &  t4,
double  jac[4][4],
double &  det,
int  imap 
)
protected

Definition at line 3257 of file ComponentFieldMap.cc.

3259 {
3260
3261 if (debug) {
3262 std::cout << m_className << "::Coordinates13:\n";
3263 std::cout << " Point (" << x << ", " << y << ", " << z << ")\n";
3264 }
3265
3266 // Failure flag
3267 int ifail = 1;
3268
3269 // Provisional values
3270 t1 = t2 = t3 = t4 = 0.;
3271
3272 // Set tolerance parameter.
3273 double f = 0.5;
3274
3275 // Make a first order approximation.
3276 int rc = Coordinates12(x, y, z, t1, t2, t3, t4, imap);
3277 if (rc > 0) {
3278 if (debug) {
3279 std::cout << m_className << "::Coordinates13:\n";
3280 std::cout << " Failure to obtain linear estimate of isoparametric "
3281 "coordinates\n";
3282 std::cout << " in element " << imap << ".\n";
3283 }
3284 return ifail;
3285 }
3286 if (t1 < -f || t2 < -f || t3 < -f || t4 < -f || t1 > 1 + f || t2 > 1 + f ||
3287 t3 > 1 + f || t4 > 1 + f) {
3288 if (debug) {
3289 std::cout << m_className << "::Coordinates13:\n";
3290 std::cout << " Linear isoparametric coordinates more than\n";
3291 std::cout << " f (" << f << ") out of range in element " << imap
3292 << ".\n";
3293 }
3294 ifail = 0;
3295 return ifail;
3296 }
3297
3298 // Start iteration.
3299 double td1 = t1, td2 = t2, td3 = t3, td4 = t4;
3300 if (debug) {
3301 std::cout << m_className << "::Coordinates13:\n";
3302 std::cout << " Iteration starts at (t1,t2,t3,t4) = (" << td1 << ", "
3303 << td2 << ", " << td3 << ", " << td4 << ").\n";
3304 }
3305 // Loop
3306 bool converged = false;
3307 double diff[4], corr[4];
3308 for (int iter = 0; iter < 10; iter++) {
3309 if (debug) {
3310 std::cout << m_className << "::Coordinates13:\n";
3311 std::cout << " Iteration " << iter << ": (t1,t2,t3,t4) = (" << td1
3312 << ", " << td2 << ", " << td3 << ", " << td4 << ").\n";
3313 }
3314 // Re-compute the (x,y,z) position for this coordinate.
3315 double xr = nodes[elements[imap].emap[0]].x * td1 * (2 * td1 - 1) +
3316 nodes[elements[imap].emap[1]].x * td2 * (2 * td2 - 1) +
3317 nodes[elements[imap].emap[2]].x * td3 * (2 * td3 - 1) +
3318 nodes[elements[imap].emap[3]].x * td4 * (2 * td4 - 1) +
3319 nodes[elements[imap].emap[4]].x * 4 * td1 * td2 +
3320 nodes[elements[imap].emap[5]].x * 4 * td1 * td3 +
3321 nodes[elements[imap].emap[6]].x * 4 * td1 * td4 +
3322 nodes[elements[imap].emap[7]].x * 4 * td2 * td3 +
3323 nodes[elements[imap].emap[8]].x * 4 * td2 * td4 +
3324 nodes[elements[imap].emap[9]].x * 4 * td3 * td4;
3325 double yr = nodes[elements[imap].emap[0]].y * td1 * (2 * td1 - 1) +
3326 nodes[elements[imap].emap[1]].y * td2 * (2 * td2 - 1) +
3327 nodes[elements[imap].emap[2]].y * td3 * (2 * td3 - 1) +
3328 nodes[elements[imap].emap[3]].y * td4 * (2 * td4 - 1) +
3329 nodes[elements[imap].emap[4]].y * 4 * td1 * td2 +
3330 nodes[elements[imap].emap[5]].y * 4 * td1 * td3 +
3331 nodes[elements[imap].emap[6]].y * 4 * td1 * td4 +
3332 nodes[elements[imap].emap[7]].y * 4 * td2 * td3 +
3333 nodes[elements[imap].emap[8]].y * 4 * td2 * td4 +
3334 nodes[elements[imap].emap[9]].y * 4 * td3 * td4;
3335 double zr = nodes[elements[imap].emap[0]].z * td1 * (2 * td1 - 1) +
3336 nodes[elements[imap].emap[1]].z * td2 * (2 * td2 - 1) +
3337 nodes[elements[imap].emap[2]].z * td3 * (2 * td3 - 1) +
3338 nodes[elements[imap].emap[3]].z * td4 * (2 * td4 - 1) +
3339 nodes[elements[imap].emap[4]].z * 4 * td1 * td2 +
3340 nodes[elements[imap].emap[5]].z * 4 * td1 * td3 +
3341 nodes[elements[imap].emap[6]].z * 4 * td1 * td4 +
3342 nodes[elements[imap].emap[7]].z * 4 * td2 * td3 +
3343 nodes[elements[imap].emap[8]].z * 4 * td2 * td4 +
3344 nodes[elements[imap].emap[9]].z * 4 * td3 * td4;
3345 double sr = td1 + td2 + td3 + td4;
3346
3347 // Compute the Jacobian.
3348 Jacobian13(imap, td1, td2, td3, td4, det, jac);
3349 // Compute the difference vector.
3350 diff[0] = 1 - sr;
3351 diff[1] = x - xr;
3352 diff[2] = y - yr;
3353 diff[3] = z - zr;
3354
3355 // Update the estimate.
3356 for (int l = 0; l < 4; ++l) {
3357 corr[l] = 0;
3358 for (int k = 0; k < 4; ++k) {
3359 corr[l] += jac[l][k] * diff[k] / det;
3360 }
3361 }
3362
3363 // Debugging
3364 if (debug) {
3365 std::cout << m_className << "::Coordinates13:\n";
3366 std::cout << " Difference vector: (1, x, y, z) = (" << diff[0]
3367 << ", " << diff[1] << ", " << diff[2] << ", " << diff[3]
3368 << ").\n";
3369 std::cout << " Correction vector: (t1,t2,t3,t4) = (" << corr[0]
3370 << ", " << corr[1] << ", " << corr[2] << ", " << corr[3]
3371 << ").\n";
3372 }
3373
3374 // Update the vector.
3375 td1 += corr[0];
3376 td2 += corr[1];
3377 td3 += corr[2];
3378 td4 += corr[3];
3379
3380 // Check for convergence.
3381 if (fabs(corr[0]) < 1.0e-5 && fabs(corr[1]) < 1.0e-5 &&
3382 fabs(corr[2]) < 1.0e-5 && fabs(corr[3]) < 1.0e-5) {
3383 if (debug) {
3384 std::cout << m_className << "::Coordinates13:\n";
3385 std::cout << " Convergence reached.\n";
3386 }
3387 converged = true;
3388 break;
3389 }
3390 }
3391
3392 // No convergence reached.
3393 if (!converged) {
3394 double xmin, ymin, zmin, xmax, ymax, zmax;
3395 xmin = nodes[elements[imap].emap[0]].x;
3396 xmax = nodes[elements[imap].emap[0]].x;
3397 if (nodes[elements[imap].emap[1]].x < xmin) {
3398 xmin = nodes[elements[imap].emap[1]].x;
3399 }
3400 if (nodes[elements[imap].emap[1]].x > xmax) {
3401 xmax = nodes[elements[imap].emap[1]].x;
3402 }
3403 if (nodes[elements[imap].emap[2]].x < xmin) {
3404 xmin = nodes[elements[imap].emap[2]].x;
3405 }
3406 if (nodes[elements[imap].emap[2]].x > xmax) {
3407 xmax = nodes[elements[imap].emap[2]].x;
3408 }
3409 if (nodes[elements[imap].emap[3]].x < xmin) {
3410 xmin = nodes[elements[imap].emap[3]].x;
3411 }
3412 if (nodes[elements[imap].emap[3]].x > xmax) {
3413 xmax = nodes[elements[imap].emap[3]].x;
3414 }
3415 ymin = nodes[elements[imap].emap[0]].y;
3416 ymax = nodes[elements[imap].emap[0]].y;
3417 if (nodes[elements[imap].emap[1]].y < ymin) {
3418 ymin = nodes[elements[imap].emap[1]].y;
3419 }
3420 if (nodes[elements[imap].emap[1]].y > ymax) {
3421 ymax = nodes[elements[imap].emap[1]].y;
3422 }
3423 if (nodes[elements[imap].emap[2]].y < ymin) {
3424 ymin = nodes[elements[imap].emap[2]].y;
3425 }
3426 if (nodes[elements[imap].emap[2]].y > ymax) {
3427 ymax = nodes[elements[imap].emap[2]].y;
3428 }
3429 if (nodes[elements[imap].emap[3]].y < ymin) {
3430 ymin = nodes[elements[imap].emap[3]].y;
3431 }
3432 if (nodes[elements[imap].emap[3]].y > ymax) {
3433 ymax = nodes[elements[imap].emap[3]].y;
3434 }
3435 zmin = nodes[elements[imap].emap[0]].z;
3436 zmax = nodes[elements[imap].emap[0]].z;
3437 if (nodes[elements[imap].emap[1]].z < zmin) {
3438 zmin = nodes[elements[imap].emap[1]].z;
3439 }
3440 if (nodes[elements[imap].emap[1]].z > zmax) {
3441 zmax = nodes[elements[imap].emap[1]].z;
3442 }
3443 if (nodes[elements[imap].emap[2]].z < zmin) {
3444 zmin = nodes[elements[imap].emap[2]].z;
3445 }
3446 if (nodes[elements[imap].emap[2]].z > zmax) {
3447 zmax = nodes[elements[imap].emap[2]].z;
3448 }
3449 if (nodes[elements[imap].emap[3]].z < zmin) {
3450 zmin = nodes[elements[imap].emap[3]].z;
3451 }
3452 if (nodes[elements[imap].emap[3]].z > zmax) {
3453 zmax = nodes[elements[imap].emap[3]].z;
3454 }
3455
3456 if (x >= xmin && x <= xmax && y >= ymin && y <= ymax && z >= zmin &&
3457 z <= zmax) {
3458 std::cout << m_className << "::Coordinates13:\n";
3459 std::cout << " No convergence achieved "
3460 << "when refining internal isoparametric coordinates\n";
3461 std::cout << " in element " << imap << " at position (" << x << ", "
3462 << y << ", " << z << ").\n";
3463 t1 = t2 = t3 = t4 = -1;
3464 return ifail;
3465 }
3466 }
3467
3468 // Convergence reached.
3469 t1 = td1;
3470 t2 = td2;
3471 t3 = td3;
3472 t4 = td4;
3473 if (debug) {
3474 std::cout << m_className << "::Coordinates13:\n";
3475 std::cout << " Convergence reached at (t1, t2, t3, t4) = (" << t1 << ", "
3476 << t2 << ", " << t3 << ", " << t4 << ").\n";
3477 }
3478
3479 // For debugging purposes, show position.
3480 if (debug) {
3481 // Re-compute the (x,y,z) position for this coordinate.
3482 double xr = nodes[elements[imap].emap[0]].x * td1 * (2 * td1 - 1) +
3483 nodes[elements[imap].emap[1]].x * td2 * (2 * td2 - 1) +
3484 nodes[elements[imap].emap[2]].x * td3 * (2 * td3 - 1) +
3485 nodes[elements[imap].emap[3]].x * td4 * (2 * td4 - 1) +
3486 nodes[elements[imap].emap[4]].x * 4 * td1 * td2 +
3487 nodes[elements[imap].emap[5]].x * 4 * td1 * td3 +
3488 nodes[elements[imap].emap[6]].x * 4 * td1 * td4 +
3489 nodes[elements[imap].emap[7]].x * 4 * td2 * td3 +
3490 nodes[elements[imap].emap[8]].x * 4 * td2 * td4 +
3491 nodes[elements[imap].emap[9]].x * 4 * td3 * td4;
3492 double yr = nodes[elements[imap].emap[0]].y * td1 * (2 * td1 - 1) +
3493 nodes[elements[imap].emap[1]].y * td2 * (2 * td2 - 1) +
3494 nodes[elements[imap].emap[2]].y * td3 * (2 * td3 - 1) +
3495 nodes[elements[imap].emap[3]].y * td4 * (2 * td4 - 1) +
3496 nodes[elements[imap].emap[4]].y * 4 * td1 * td2 +
3497 nodes[elements[imap].emap[5]].y * 4 * td1 * td3 +
3498 nodes[elements[imap].emap[6]].y * 4 * td1 * td4 +
3499 nodes[elements[imap].emap[7]].y * 4 * td2 * td3 +
3500 nodes[elements[imap].emap[8]].y * 4 * td2 * td4 +
3501 nodes[elements[imap].emap[9]].y * 4 * td3 * td4;
3502 double zr = nodes[elements[imap].emap[0]].z * td1 * (2 * td1 - 1) +
3503 nodes[elements[imap].emap[1]].z * td2 * (2 * td2 - 1) +
3504 nodes[elements[imap].emap[2]].z * td3 * (2 * td3 - 1) +
3505 nodes[elements[imap].emap[3]].z * td4 * (2 * td4 - 1) +
3506 nodes[elements[imap].emap[4]].z * 4 * td1 * td2 +
3507 nodes[elements[imap].emap[5]].z * 4 * td1 * td3 +
3508 nodes[elements[imap].emap[6]].z * 4 * td1 * td4 +
3509 nodes[elements[imap].emap[7]].z * 4 * td2 * td3 +
3510 nodes[elements[imap].emap[8]].z * 4 * td2 * td4 +
3511 nodes[elements[imap].emap[9]].z * 4 * td3 * td4;
3512 double sr = td1 + td2 + td3 + td4;
3513 std::cout << m_className << "::Coordinates13:\n";
3514 std::cout << " Position requested: (" << x << ", " << y << ", " << z
3515 << ")\n";
3516 std::cout << " Reconstructed: (" << xr << ", " << yr << ", "
3517 << zr << ")\n";
3518 std::cout << " Difference: (" << x - xr << ", " << y - yr
3519 << ", " << z - zr << ")\n";
3520 std::cout << " Checksum - 1: " << sr - 1 << "\n";
3521 }
3522
3523 // Success
3524 ifail = 0;
3525 return ifail;
3526}
fabs
DoubleAc fabs(const DoubleAc &f)
Definition: DoubleAc.h:616
Garfield::ComponentFieldMap::Coordinates12
int Coordinates12(double x, double y, double z, double &t1, double &t2, double &t3, double &t4, int imap)
Definition: ComponentFieldMap.cc:3005
Garfield::ComponentFieldMap::Jacobian13
void Jacobian13(int i, double t, double u, double v, double w, double &det, double jac[4][4])
Definition: ComponentFieldMap.cc:1233

Referenced by FindElement13().

◆ Coordinates3()

int Garfield::ComponentFieldMap::Coordinates3 ( double  x,
double  y,
double  z,
double &  t1,
double &  t2,
double &  t3,
double &  t4,
double  jac[4][4],
double &  det,
int  imap 
)
protected

Definition at line 2226 of file ComponentFieldMap.cc.

2228 {
2229
2230 if (debug) {
2231 std::cout << m_className << "::Coordinates3:\n";
2232 std::cout << " Point (" << x << ", " << y << ", " << z << ")\n";
2233 }
2234
2235 // Failure flag
2236 int ifail = 1;
2237
2238 // Provisional values
2239 t1 = t2 = t3 = t4 = 0;
2240
2241 // Make a first order approximation, using the linear triangle.
2242 double tt1 =
2243 (x - nodes[elements[imap].emap[1]].x) *
2244 (nodes[elements[imap].emap[2]].y - nodes[elements[imap].emap[1]].y) -
2245 (y - nodes[elements[imap].emap[1]].y) *
2246 (nodes[elements[imap].emap[2]].x - nodes[elements[imap].emap[1]].x);
2247 double tt2 =
2248 (x - nodes[elements[imap].emap[2]].x) *
2249 (nodes[elements[imap].emap[0]].y - nodes[elements[imap].emap[2]].y) -
2250 (y - nodes[elements[imap].emap[2]].y) *
2251 (nodes[elements[imap].emap[0]].x - nodes[elements[imap].emap[2]].x);
2252 double tt3 =
2253 (x - nodes[elements[imap].emap[0]].x) *
2254 (nodes[elements[imap].emap[1]].y - nodes[elements[imap].emap[0]].y) -
2255 (y - nodes[elements[imap].emap[0]].y) *
2256 (nodes[elements[imap].emap[1]].x - nodes[elements[imap].emap[0]].x);
2257 if ((nodes[elements[imap].emap[0]].x - nodes[elements[imap].emap[1]].x) *
2258 (nodes[elements[imap].emap[2]].y -
2259 nodes[elements[imap].emap[1]].y) -
2260 (nodes[elements[imap].emap[2]].x -
2261 nodes[elements[imap].emap[1]].x) *
2262 (nodes[elements[imap].emap[0]].y -
2263 nodes[elements[imap].emap[1]].y) ==
2264 0 ||
2265 (nodes[elements[imap].emap[1]].x - nodes[elements[imap].emap[2]].x) *
2266 (nodes[elements[imap].emap[0]].y -
2267 nodes[elements[imap].emap[2]].y) -
2268 (nodes[elements[imap].emap[0]].x -
2269 nodes[elements[imap].emap[2]].x) *
2270 (nodes[elements[imap].emap[1]].y -
2271 nodes[elements[imap].emap[2]].y) ==
2272 0 ||
2273 (nodes[elements[imap].emap[2]].x - nodes[elements[imap].emap[0]].x) *
2274 (nodes[elements[imap].emap[1]].y -
2275 nodes[elements[imap].emap[0]].y) -
2276 (nodes[elements[imap].emap[1]].x -
2277 nodes[elements[imap].emap[0]].x) *
2278 (nodes[elements[imap].emap[2]].y -
2279 nodes[elements[imap].emap[0]].y) ==
2280 0) {
2281 std::cerr << m_className << "::Coordinates3:\n";
2282 std::cerr << " Calculation of linear coordinates failed; abandoned.\n";
2283 return ifail;
2284 } else {
2285 t1 = tt1 /
2286 ((nodes[elements[imap].emap[0]].x - nodes[elements[imap].emap[1]].x) *
2287 (nodes[elements[imap].emap[2]].y -
2288 nodes[elements[imap].emap[1]].y) -
2289 (nodes[elements[imap].emap[2]].x - nodes[elements[imap].emap[1]].x) *
2290 (nodes[elements[imap].emap[0]].y -
2291 nodes[elements[imap].emap[1]].y));
2292 t2 = tt2 /
2293 ((nodes[elements[imap].emap[1]].x - nodes[elements[imap].emap[2]].x) *
2294 (nodes[elements[imap].emap[0]].y -
2295 nodes[elements[imap].emap[2]].y) -
2296 (nodes[elements[imap].emap[0]].x - nodes[elements[imap].emap[2]].x) *
2297 (nodes[elements[imap].emap[1]].y -
2298 nodes[elements[imap].emap[2]].y));
2299 t3 = tt3 /
2300 ((nodes[elements[imap].emap[2]].x - nodes[elements[imap].emap[0]].x) *
2301 (nodes[elements[imap].emap[1]].y -
2302 nodes[elements[imap].emap[0]].y) -
2303 (nodes[elements[imap].emap[1]].x - nodes[elements[imap].emap[0]].x) *
2304 (nodes[elements[imap].emap[2]].y -
2305 nodes[elements[imap].emap[0]].y));
2306 t4 = 0;
2307 }
2308
2309 // Start iterative refinement
2310 double td1 = t1, td2 = t2, td3 = t3;
2311 if (debug) {
2312 std::cout << m_className << "::Coordinates3:\n";
2313 std::cout << " Linear estimate: (u, v, w) = (" << td1 << ", " << td2
2314 << ", " << td3 << "), sum = " << td1 + td2 + td3 << ".\n";
2315 }
2316
2317 // Loop
2318 bool converged = false;
2319 double diff[3] = {0., 0., 0.};
2320 double corr[3] = {0., 0., 0.};
2321 for (int iter = 0; iter < 10; iter++) {
2322 if (debug) {
2323 std::cout << m_className << "::Coordinates3:\n";
2324 std::cout << " Iteration " << iter << ": (u, v, w) = (" << td1
2325 << ", " << td2 << ", " << td3 << "), sum = " << td1 + td2 + td3
2326 << "\n";
2327 }
2328 // Re-compute the (x,y,z) position for this coordinate.
2329 double xr = nodes[elements[imap].emap[0]].x * td1 * (2 * td1 - 1) +
2330 nodes[elements[imap].emap[1]].x * td2 * (2 * td2 - 1) +
2331 nodes[elements[imap].emap[2]].x * td3 * (2 * td3 - 1) +
2332 nodes[elements[imap].emap[3]].x * 4 * td1 * td2 +
2333 nodes[elements[imap].emap[4]].x * 4 * td1 * td3 +
2334 nodes[elements[imap].emap[5]].x * 4 * td2 * td3;
2335 double yr = nodes[elements[imap].emap[0]].y * td1 * (2 * td1 - 1) +
2336 nodes[elements[imap].emap[1]].y * td2 * (2 * td2 - 1) +
2337 nodes[elements[imap].emap[2]].y * td3 * (2 * td3 - 1) +
2338 nodes[elements[imap].emap[3]].y * 4 * td1 * td2 +
2339 nodes[elements[imap].emap[4]].y * 4 * td1 * td3 +
2340 nodes[elements[imap].emap[5]].y * 4 * td2 * td3;
2341 double sr = td1 + td2 + td3;
2342 // Compute the Jacobian
2343 Jacobian3(imap, td1, td2, td3, det, jac);
2344 // Compute the difference vector
2345 diff[0] = 1 - sr;
2346 diff[1] = x - xr;
2347 diff[2] = y - yr;
2348 // Update the estimate
2349 for (int l = 0; l < 3; l++) {
2350 corr[l] = 0;
2351 for (int k = 0; k < 3; k++) {
2352 corr[l] += jac[l][k] * diff[k] / det;
2353 }
2354 }
2355 // Debugging
2356 if (debug) {
2357 std::cout << m_className << "::Coordinates3:\n";
2358 std::cout << " Difference vector: (1, x, y) = (" << diff[0] << ", "
2359 << diff[1] << ", " << diff[2] << ").\n";
2360 std::cout << " Correction vector: (u, v, w) = (" << corr[0] << ", "
2361 << corr[1] << ", " << corr[2] << ").\n";
2362 }
2363 // Update the vector.
2364 td1 += corr[0];
2365 td2 += corr[1];
2366 td3 += corr[2];
2367 // Check for convergence.
2368 if (fabs(corr[0]) < 1.0e-5 && fabs(corr[1]) < 1.0e-5 &&
2369 fabs(corr[2]) < 1.0e-5) {
2370 if (debug) {
2371 std::cout << m_className << "::Coordinates3:\n";
2372 std::cout << " Convergence reached.\n";
2373 }
2374 converged = true;
2375 break;
2376 }
2377 }
2378 // No convergence reached
2379 if (!converged) {
2380 double xmin, ymin, xmax, ymax;
2381 xmin = nodes[elements[imap].emap[0]].x;
2382 xmax = nodes[elements[imap].emap[0]].x;
2383 if (nodes[elements[imap].emap[1]].x < xmin) {
2384 xmin = nodes[elements[imap].emap[1]].x;
2385 }
2386 if (nodes[elements[imap].emap[1]].x > xmax) {
2387 xmax = nodes[elements[imap].emap[1]].x;
2388 }
2389 if (nodes[elements[imap].emap[2]].x < xmin) {
2390 xmin = nodes[elements[imap].emap[2]].x;
2391 }
2392 if (nodes[elements[imap].emap[2]].x > xmax) {
2393 xmax = nodes[elements[imap].emap[2]].x;
2394 }
2395 ymin = nodes[elements[imap].emap[0]].y;
2396 ymax = nodes[elements[imap].emap[0]].y;
2397 if (nodes[elements[imap].emap[1]].y < ymin) {
2398 ymin = nodes[elements[imap].emap[1]].y;
2399 }
2400 if (nodes[elements[imap].emap[1]].y > ymax) {
2401 ymax = nodes[elements[imap].emap[1]].y;
2402 }
2403 if (nodes[elements[imap].emap[2]].y < ymin) {
2404 ymin = nodes[elements[imap].emap[2]].y;
2405 }
2406 if (nodes[elements[imap].emap[2]].y > ymax) {
2407 ymax = nodes[elements[imap].emap[2]].y;
2408 }
2409
2410 if (x >= xmin && x <= xmax && y >= ymin && y <= ymax) {
2411 std::cout << m_className << "::Coordinates3:\n";
2412 std::cout << " No convergence achieved "
2413 << "when refining internal isoparametric coordinates\n";
2414 std::cout << " in element " << imap << " at position (" << x << ", "
2415 << y << ").\n";
2416 t1 = t2 = t3 = t4 = 0;
2417 return ifail;
2418 }
2419 }
2420
2421 // Convergence reached.
2422 t1 = td1;
2423 t2 = td2;
2424 t3 = td3;
2425 t4 = 0;
2426 if (debug) {
2427 std::cout << m_className << "::Coordinates3:\n";
2428 std::cout << " Convergence reached at (t1, t2, t3) = (" << t1 << ", "
2429 << t2 << ", " << t3 << ").\n";
2430 }
2431
2432 // For debugging purposes, show position
2433 if (debug) {
2434 double xr = nodes[elements[imap].emap[0]].x * td1 * (2 * td1 - 1) +
2435 nodes[elements[imap].emap[1]].x * td2 * (2 * td2 - 1) +
2436 nodes[elements[imap].emap[2]].x * td3 * (2 * td3 - 1) +
2437 nodes[elements[imap].emap[3]].x * 4 * td1 * td2 +
2438 nodes[elements[imap].emap[4]].x * 4 * td1 * td3 +
2439 nodes[elements[imap].emap[5]].x * 4 * td2 * td3;
2440 double yr = nodes[elements[imap].emap[0]].y * td1 * (2 * td1 - 1) +
2441 nodes[elements[imap].emap[1]].y * td2 * (2 * td2 - 1) +
2442 nodes[elements[imap].emap[2]].y * td3 * (2 * td3 - 1) +
2443 nodes[elements[imap].emap[3]].y * 4 * td1 * td2 +
2444 nodes[elements[imap].emap[4]].y * 4 * td1 * td3 +
2445 nodes[elements[imap].emap[5]].y * 4 * td2 * td3;
2446 double sr = td1 + td2 + td3;
2447 std::cout << m_className << "::Coordinates3:\n";
2448 std::cout << " Position requested: (" << x << ", " << y << ")\n";
2449 std::cout << " Reconstructed: (" << xr << ", " << yr << ")\n";
2450 std::cout << " Difference: (" << x - xr << ", " << y - yr
2451 << ")\n";
2452 std::cout << " Checksum - 1: " << sr - 1 << "\n";
2453 }
2454
2455 // Success
2456 ifail = 0;
2457 return ifail;
2458}
Garfield::ComponentFieldMap::Jacobian3
void Jacobian3(int i, double u, double v, double w, double &det, double jac[4][4])
Definition: ComponentFieldMap.cc:558

Referenced by FindElement5().

◆ Coordinates4()

int Garfield::ComponentFieldMap::Coordinates4 ( double  x,
double  y,
double  z,
double &  t1,
double &  t2,
double &  t3,
double &  t4,
double  jac[4][4],
double &  det,
int  imap 
)
protected

Definition at line 2460 of file ComponentFieldMap.cc.

2462 {
2463
2464 // Debugging
2465 if (debug) {
2466 std::cout << m_className << "::Coordinates4:\n";
2467 std::cout << " Point (" << x << ", " << y << ", " << z << ")\n";
2468 }
2469
2470 // Failure flag
2471 int ifail = 1;
2472
2473 // Provisional values
2474 t1 = t2 = t3 = t4 = 0.;
2475
2476 // Compute determinant.
2477 det =
2478 -(-((nodes[elements[imap].emap[0]].x - nodes[elements[imap].emap[3]].x) *
2479 (nodes[elements[imap].emap[1]].y - nodes[elements[imap].emap[2]].y)) +
2480 (nodes[elements[imap].emap[1]].x - nodes[elements[imap].emap[2]].x) *
2481 (nodes[elements[imap].emap[0]].y -
2482 nodes[elements[imap].emap[3]].y)) *
2483 (2 * x * (-nodes[elements[imap].emap[0]].y +
2484 nodes[elements[imap].emap[1]].y +
2485 nodes[elements[imap].emap[2]].y -
2486 nodes[elements[imap].emap[3]].y) -
2487 (nodes[elements[imap].emap[0]].x + nodes[elements[imap].emap[3]].x) *
2488 (nodes[elements[imap].emap[1]].y +
2489 nodes[elements[imap].emap[2]].y - 2 * y) +
2490 nodes[elements[imap].emap[1]].x *
2491 (nodes[elements[imap].emap[0]].y +
2492 nodes[elements[imap].emap[3]].y - 2 * y) +
2493 nodes[elements[imap].emap[2]].x *
2494 (nodes[elements[imap].emap[0]].y +
2495 nodes[elements[imap].emap[3]].y - 2 * y)) +
2496 pow(-(nodes[elements[imap].emap[0]].x * nodes[elements[imap].emap[1]].y) +
2497 nodes[elements[imap].emap[3]].x *
2498 nodes[elements[imap].emap[2]].y -
2499 nodes[elements[imap].emap[2]].x *
2500 nodes[elements[imap].emap[3]].y +
2501 x * (-nodes[elements[imap].emap[0]].y +
2502 nodes[elements[imap].emap[1]].y -
2503 nodes[elements[imap].emap[2]].y +
2504 nodes[elements[imap].emap[3]].y) +
2505 nodes[elements[imap].emap[1]].x *
2506 (nodes[elements[imap].emap[0]].y - y) +
2507 (nodes[elements[imap].emap[0]].x +
2508 nodes[elements[imap].emap[2]].x -
2509 nodes[elements[imap].emap[3]].x) *
2510 y,
2511 2);
2512
2513 // Check that the determinant is non-negative
2514 // (this can happen if the point is out of range).
2515 if (det < 0) {
2516 if (debug) {
2517 std::cerr << m_className << "::Coordinates4:\n";
2518 std::cerr << " No solution found for isoparametric coordinates\n";
2519 std::cerr << " in element " << imap << " because the determinant "
2520 << det << " is < 0.\n";
2521 }
2522 return ifail;
2523 }
2524
2525 // Vector products for evaluation of T1.
2526 double prod =
2527 ((nodes[elements[imap].emap[2]].x - nodes[elements[imap].emap[3]].x) *
2528 (nodes[elements[imap].emap[0]].y - nodes[elements[imap].emap[1]].y) -
2529 (nodes[elements[imap].emap[0]].x - nodes[elements[imap].emap[1]].x) *
2530 (nodes[elements[imap].emap[2]].y - nodes[elements[imap].emap[3]].y));
2531 if (prod * prod >
2532 1.0e-12 *
2533 ((nodes[elements[imap].emap[0]].x - nodes[elements[imap].emap[1]].x) *
2534 (nodes[elements[imap].emap[0]].x -
2535 nodes[elements[imap].emap[1]].x) +
2536 (nodes[elements[imap].emap[0]].y - nodes[elements[imap].emap[1]].y) *
2537 (nodes[elements[imap].emap[0]].y -
2538 nodes[elements[imap].emap[1]].y)) *
2539 ((nodes[elements[imap].emap[2]].x - nodes[elements[imap].emap[3]].x) *
2540 (nodes[elements[imap].emap[2]].x -
2541 nodes[elements[imap].emap[3]].x) +
2542 (nodes[elements[imap].emap[2]].y - nodes[elements[imap].emap[3]].y) *
2543 (nodes[elements[imap].emap[2]].y -
2544 nodes[elements[imap].emap[3]].y))) {
2545 t1 = (-(nodes[elements[imap].emap[3]].x * nodes[elements[imap].emap[0]].y) +
2546 x * nodes[elements[imap].emap[0]].y +
2547 nodes[elements[imap].emap[2]].x * nodes[elements[imap].emap[1]].y -
2548 x * nodes[elements[imap].emap[1]].y -
2549 nodes[elements[imap].emap[1]].x * nodes[elements[imap].emap[2]].y +
2550 x * nodes[elements[imap].emap[2]].y +
2551 nodes[elements[imap].emap[0]].x * nodes[elements[imap].emap[3]].y -
2552 x * nodes[elements[imap].emap[3]].y -
2553 nodes[elements[imap].emap[0]].x * y +
2554 nodes[elements[imap].emap[1]].x * y -
2555 nodes[elements[imap].emap[2]].x * y +
2556 nodes[elements[imap].emap[3]].x * y + sqrt(det)) /
2557 prod;
2558 } else {
2559 double xp =
2560 nodes[elements[imap].emap[0]].y - nodes[elements[imap].emap[1]].y;
2561 double yp =
2562 nodes[elements[imap].emap[1]].x - nodes[elements[imap].emap[0]].x;
2563 double dn = sqrt(xp * xp + yp * yp);
2564 if (dn <= 0) {
2565 std::cerr << m_className << "::Coordinates4:\n";
2566 std::cerr << " Element " << imap
2567 << " appears to be degenerate in the 1 - 2 axis.\n";
2568 return ifail;
2569 }
2570 xp = xp / dn;
2571 yp = yp / dn;
2572 double dpoint = xp * (x - nodes[elements[imap].emap[0]].x) +
2573 yp * (y - nodes[elements[imap].emap[0]].y);
2574 double dbox = xp * (nodes[elements[imap].emap[3]].x -
2575 nodes[elements[imap].emap[0]].x) +
2576 yp * (nodes[elements[imap].emap[3]].y -
2577 nodes[elements[imap].emap[0]].y);
2578 if (dbox == 0) {
2579 std::cerr << m_className << "::Coordinates4:\n";
2580 std::cerr << " Element " << imap
2581 << " appears to be degenerate in the 1 - 3 axis.\n";
2582 return ifail;
2583 }
2584 double t = -1 + 2 * dpoint / dbox;
2585 double xt1 = nodes[elements[imap].emap[0]].x +
2586 0.5 * (t + 1) * (nodes[elements[imap].emap[3]].x -
2587 nodes[elements[imap].emap[0]].x);
2588 double yt1 = nodes[elements[imap].emap[0]].y +
2589 0.5 * (t + 1) * (nodes[elements[imap].emap[3]].y -
2590 nodes[elements[imap].emap[0]].y);
2591 double xt2 = nodes[elements[imap].emap[1]].x +
2592 0.5 * (t + 1) * (nodes[elements[imap].emap[2]].x -
2593 nodes[elements[imap].emap[1]].x);
2594 double yt2 = nodes[elements[imap].emap[1]].y +
2595 0.5 * (t + 1) * (nodes[elements[imap].emap[2]].y -
2596 nodes[elements[imap].emap[1]].y);
2597 dn = (xt1 - xt2) * (xt1 - xt2) + (yt1 - yt2) * (yt1 - yt2);
2598 if (dn <= 0) {
2599 std::cout << m_className << "::Coordinates4:\n";
2600 std::cout << " Coordinate requested at convergence point of element "
2601 << imap << ".\n";
2602 return ifail;
2603 }
2604 t1 = -1 + 2 * ((x - xt1) * (xt2 - xt1) + (y - yt1) * (yt2 - yt1)) / dn;
2605 }
2606
2607 // Vector products for evaluation of T2.
2608 prod =
2609 ((nodes[elements[imap].emap[0]].x - nodes[elements[imap].emap[3]].x) *
2610 (nodes[elements[imap].emap[1]].y - nodes[elements[imap].emap[2]].y) -
2611 (nodes[elements[imap].emap[1]].x - nodes[elements[imap].emap[2]].x) *
2612 (nodes[elements[imap].emap[0]].y - nodes[elements[imap].emap[3]].y));
2613 if (prod * prod >
2614 1.0e-12 *
2615 ((nodes[elements[imap].emap[0]].x - nodes[elements[imap].emap[3]].x) *
2616 (nodes[elements[imap].emap[0]].x -
2617 nodes[elements[imap].emap[3]].x) +
2618 (nodes[elements[imap].emap[0]].y - nodes[elements[imap].emap[3]].y) *
2619 (nodes[elements[imap].emap[0]].y -
2620 nodes[elements[imap].emap[3]].y)) *
2621 ((nodes[elements[imap].emap[1]].x - nodes[elements[imap].emap[2]].x) *
2622 (nodes[elements[imap].emap[1]].x -
2623 nodes[elements[imap].emap[2]].x) +
2624 (nodes[elements[imap].emap[1]].y - nodes[elements[imap].emap[2]].y) *
2625 (nodes[elements[imap].emap[1]].y -
2626 nodes[elements[imap].emap[2]].y))) {
2627 t2 = (-(nodes[elements[imap].emap[1]].x * nodes[elements[imap].emap[0]].y) +
2628 x * nodes[elements[imap].emap[0]].y +
2629 nodes[elements[imap].emap[0]].x * nodes[elements[imap].emap[1]].y -
2630 x * nodes[elements[imap].emap[1]].y -
2631 nodes[elements[imap].emap[3]].x * nodes[elements[imap].emap[2]].y +
2632 x * nodes[elements[imap].emap[2]].y +
2633 nodes[elements[imap].emap[2]].x * nodes[elements[imap].emap[3]].y -
2634 x * nodes[elements[imap].emap[3]].y -
2635 nodes[elements[imap].emap[0]].x * y +
2636 nodes[elements[imap].emap[1]].x * y -
2637 nodes[elements[imap].emap[2]].x * y +
2638 nodes[elements[imap].emap[3]].x * y - sqrt(det)) /
2639 prod;
2640 } else {
2641 double xp =
2642 nodes[elements[imap].emap[0]].y - nodes[elements[imap].emap[3]].y;
2643 double yp =
2644 nodes[elements[imap].emap[3]].x - nodes[elements[imap].emap[0]].x;
2645 double dn = sqrt(xp * xp + yp * yp);
2646 if (dn <= 0) {
2647 std::cerr << m_className << "Coordinates4:\n";
2648 std::cerr << " Element " << imap
2649 << " appears to be degenerate in the 1 - 4 axis.\n";
2650 return ifail;
2651 }
2652 xp = xp / dn;
2653 yp = yp / dn;
2654 double dpoint = xp * (x - nodes[elements[imap].emap[0]].x) +
2655 yp * (y - nodes[elements[imap].emap[0]].y);
2656 double dbox = xp * (nodes[elements[imap].emap[1]].x -
2657 nodes[elements[imap].emap[0]].x) +
2658 yp * (nodes[elements[imap].emap[1]].y -
2659 nodes[elements[imap].emap[0]].y);
2660 if (dbox == 0) {
2661 std::cerr << m_className << "::Coordinates4:\n";
2662 std::cerr << " Element " << imap
2663 << " appears to be degenerate in the 1 - 2 axis.\n";
2664 return ifail;
2665 }
2666 double t = -1 + 2 * dpoint / dbox;
2667 double xt1 = nodes[elements[imap].emap[0]].x +
2668 0.5 * (t + 1) * (nodes[elements[imap].emap[1]].x -
2669 nodes[elements[imap].emap[0]].x);
2670 double yt1 = nodes[elements[imap].emap[0]].y +
2671 0.5 * (t + 1) * (nodes[elements[imap].emap[1]].y -
2672 nodes[elements[imap].emap[0]].y);
2673 double xt2 = nodes[elements[imap].emap[3]].x +
2674 0.5 * (t + 1) * (nodes[elements[imap].emap[2]].x -
2675 nodes[elements[imap].emap[3]].x);
2676 double yt2 = nodes[elements[imap].emap[3]].y +
2677 0.5 * (t + 1) * (nodes[elements[imap].emap[2]].y -
2678 nodes[elements[imap].emap[3]].y);
2679 dn = (xt1 - xt2) * (xt1 - xt2) + (yt1 - yt2) * (yt1 - yt2);
2680 if (dn <= 0) {
2681 std::cout << m_className << "::Coordinates4:\n";
2682 std::cout << " Coordinate requested at convergence point of element "
2683 << imap << ".\n";
2684 return ifail;
2685 }
2686 t2 = -1 + 2 * ((x - xt1) * (xt2 - xt1) + (y - yt1) * (yt2 - yt1)) / dn;
2687 }
2688 if (debug) {
2689 std::cout << m_className << "::Coordinates4:\n";
2690 std::cout << " Isoparametric (u, v): (" << t1 << ", " << t2 << ").\n";
2691 }
2692
2693 // Re-compute the (x,y,z) position for this coordinate.
2694 if (debug) {
2695 double xr = nodes[elements[imap].emap[0]].x * (1 - t1) * (1 - t2) / 4 +
2696 nodes[elements[imap].emap[1]].x * (1 + t1) * (1 - t2) / 4 +
2697 nodes[elements[imap].emap[2]].x * (1 + t1) * (1 + t2) / 4 +
2698 nodes[elements[imap].emap[3]].x * (1 - t1) * (1 + t2) / 4;
2699 double yr = nodes[elements[imap].emap[0]].y * (1 - t1) * (1 - t2) / 4 +
2700 nodes[elements[imap].emap[1]].y * (1 + t1) * (1 - t2) / 4 +
2701 nodes[elements[imap].emap[2]].y * (1 + t1) * (1 + t2) / 4 +
2702 nodes[elements[imap].emap[3]].y * (1 - t1) * (1 + t2) / 4;
2703 std::cout << m_className << "::Coordinates4: \n";
2704 std::cout << " Position requested: (" << x << ", " << y << ")\n";
2705 std::cout << " Reconstructed: (" << xr << ", " << yr << ")\n";
2706 std::cout << " Difference: (" << x - xr << ", " << y - yr
2707 << ")\n";
2708 }
2709
2710 // This should have worked if we get this far.
2711 ifail = 0;
2712 return ifail;
2713 // Variable jac is not used.
2714 // The following lines are just for quieting the compiler.
2715 jac[0][0] = jac[0][1] = jac[0][2] = jac[0][3] = 0.;
2716 jac[1][0] = jac[1][1] = jac[1][2] = jac[1][3] = 0.;
2717 jac[2][0] = jac[2][1] = jac[2][2] = jac[2][3] = 0.;
2718 jac[3][0] = jac[3][1] = jac[3][2] = jac[3][3] = 0.;
2719}
pow
DoubleAc pow(const DoubleAc &f, double p)
Definition: DoubleAc.cpp:336
sqrt
DoubleAc sqrt(const DoubleAc &f)
Definition: DoubleAc.cpp:313

Referenced by Coordinates5().

◆ Coordinates5()

int Garfield::ComponentFieldMap::Coordinates5 ( double  x,
double  y,
double  z,
double &  t1,
double &  t2,
double &  t3,
double &  t4,
double  jac[4][4],
double &  det,
int  imap 
)
protected

Definition at line 2721 of file ComponentFieldMap.cc.

2723 {
2724
2725 // Debugging
2726 if (debug) {
2727 std::cout << m_className << "::Coordinates5:\n";
2728 std::cout << " Point (" << x << ", " << y << ", " << z << ")\n";
2729 }
2730
2731 // Failure flag
2732 int ifail = 1;
2733
2734 // Provisional values
2735 t1 = t2 = t3 = t4 = 0;
2736
2737 // Degenerate elements should have been treated as triangles.
2738 if (elements[imap].degenerate) {
2739 std::cerr << m_className << "::Coordinates5:\n";
2740 std::cerr << " Received degenerate element " << imap << ".\n";
2741 return ifail;
2742 }
2743
2744 // Set tolerance parameter.
2745 double f = 0.5;
2746
2747 // Make a first order approximation.
2748 int rc = Coordinates4(x, y, z, t1, t2, t3, t4, jac, det, imap);
2749 if (rc > 0) {
2750 if (debug) {
2751 std::cout << m_className << "::Coordinates5:\n";
2752 std::cout << " Failure to obtain linear estimate of isoparametric "
2753 "coordinates\n";
2754 std::cout << " in element " << imap << ".\n";
2755 }
2756 return ifail;
2757 }
2758
2759 // Check whether the point is far outside.
2760 if (t1 < -(1 + f) || t1 > (1 + f) || t2 < -(1 + f) || t2 > (1 + f)) {
2761 if (debug) {
2762 std::cout << m_className << "::Coordinates5:\n";
2763 std::cout << " Point far outside, (t1,t2) = (" << t1 << ", " << t2
2764 << ").\n";
2765 }
2766 return ifail;
2767 }
2768
2769 // Start iteration
2770 double td1 = t1, td2 = t2;
2771 if (debug) {
2772 std::cout << m_className << "::Coordinates5:\n";
2773 std::cout << " Iteration starts at (t1,t2) = (" << td1 << ", " << td2
2774 << ").\n";
2775 }
2776 // Loop
2777 bool converged = false;
2778 double diff[2], corr[2];
2779 for (int iter = 0; iter < 10; iter++) {
2780 if (debug) {
2781 std::cout << m_className << "::Coordinates5:\n";
2782 std::cout << " Iteration " << iter << ": (t1, t2) = (" << td1
2783 << ", " << td2 << ").\n";
2784 }
2785 // Re-compute the (x,y,z) position for this coordinate.
2786 double xr =
2787 nodes[elements[imap].emap[0]].x *
2788 (-(1 - td1) * (1 - td2) * (1 + td1 + td2)) / 4 +
2789 nodes[elements[imap].emap[1]].x *
2790 (-(1 + td1) * (1 - td2) * (1 - td1 + td2)) / 4 +
2791 nodes[elements[imap].emap[2]].x *
2792 (-(1 + td1) * (1 + td2) * (1 - td1 - td2)) / 4 +
2793 nodes[elements[imap].emap[3]].x *
2794 (-(1 - td1) * (1 + td2) * (1 + td1 - td2)) / 4 +
2795 nodes[elements[imap].emap[4]].x * (1 - td1) * (1 + td1) * (1 - td2) /
2796 2 +
2797 nodes[elements[imap].emap[5]].x * (1 + td1) * (1 + td2) * (1 - td2) /
2798 2 +
2799 nodes[elements[imap].emap[6]].x * (1 - td1) * (1 + td1) * (1 + td2) /
2800 2 +
2801 nodes[elements[imap].emap[7]].x * (1 - td1) * (1 + td2) * (1 - td2) / 2;
2802 double yr =
2803 nodes[elements[imap].emap[0]].y *
2804 (-(1 - td1) * (1 - td2) * (1 + td1 + td2)) / 4 +
2805 nodes[elements[imap].emap[1]].y *
2806 (-(1 + td1) * (1 - td2) * (1 - td1 + td2)) / 4 +
2807 nodes[elements[imap].emap[2]].y *
2808 (-(1 + td1) * (1 + td2) * (1 - td1 - td2)) / 4 +
2809 nodes[elements[imap].emap[3]].y *
2810 (-(1 - td1) * (1 + td2) * (1 + td1 - td2)) / 4 +
2811 nodes[elements[imap].emap[4]].y * (1 - td1) * (1 + td1) * (1 - td2) /
2812 2 +
2813 nodes[elements[imap].emap[5]].y * (1 + td1) * (1 + td2) * (1 - td2) /
2814 2 +
2815 nodes[elements[imap].emap[6]].y * (1 - td1) * (1 + td1) * (1 + td2) /
2816 2 +
2817 nodes[elements[imap].emap[7]].y * (1 - td1) * (1 + td2) * (1 - td2) / 2;
2818 // Compute the Jacobian.
2819 Jacobian5(imap, td1, td2, det, jac);
2820 // Compute the difference vector.
2821 diff[0] = x - xr;
2822 diff[1] = y - yr;
2823 // Update the estimate.
2824 for (int l = 0; l < 2; ++l) {
2825 corr[l] = 0;
2826 for (int k = 0; k < 2; ++k) {
2827 corr[l] += jac[l][k] * diff[k] / det;
2828 }
2829 }
2830 // Debugging
2831 if (debug) {
2832 std::cout << m_className << "::Coordinates5:\n";
2833 std::cout << " Difference vector: (x, y) = (" << diff[0] << ", "
2834 << diff[1] << ").\n";
2835 std::cout << " Correction vector: (t1, t2) = (" << corr[0] << ", "
2836 << corr[1] << ").\n";
2837 }
2838 // Update the vector.
2839 td1 += corr[0];
2840 td2 += corr[1];
2841 // Check for convergence.
2842 if (fabs(corr[0]) < 1.0e-5 && fabs(corr[1]) < 1.0e-5) {
2843 if (debug) {
2844 std::cout << m_className << "::Coordinates5:\n";
2845 std::cout << " Convergence reached.\n";
2846 }
2847 converged = true;
2848 break;
2849 }
2850 }
2851 // No convergence reached.
2852 if (!converged) {
2853 double xmin, ymin, xmax, ymax;
2854 xmin = nodes[elements[imap].emap[0]].x;
2855 xmax = nodes[elements[imap].emap[0]].x;
2856 if (nodes[elements[imap].emap[1]].x < xmin) {
2857 xmin = nodes[elements[imap].emap[1]].x;
2858 }
2859 if (nodes[elements[imap].emap[1]].x > xmax) {
2860 xmax = nodes[elements[imap].emap[1]].x;
2861 }
2862 if (nodes[elements[imap].emap[2]].x < xmin) {
2863 xmin = nodes[elements[imap].emap[2]].x;
2864 }
2865 if (nodes[elements[imap].emap[2]].x > xmax) {
2866 xmax = nodes[elements[imap].emap[2]].x;
2867 }
2868 if (nodes[elements[imap].emap[3]].x < xmin) {
2869 xmin = nodes[elements[imap].emap[3]].x;
2870 }
2871 if (nodes[elements[imap].emap[3]].x > xmax) {
2872 xmax = nodes[elements[imap].emap[3]].x;
2873 }
2874 if (nodes[elements[imap].emap[4]].x < xmin) {
2875 xmin = nodes[elements[imap].emap[4]].x;
2876 }
2877 if (nodes[elements[imap].emap[4]].x > xmax) {
2878 xmax = nodes[elements[imap].emap[4]].x;
2879 }
2880 if (nodes[elements[imap].emap[5]].x < xmin) {
2881 xmin = nodes[elements[imap].emap[5]].x;
2882 }
2883 if (nodes[elements[imap].emap[5]].x > xmax) {
2884 xmax = nodes[elements[imap].emap[5]].x;
2885 }
2886 if (nodes[elements[imap].emap[6]].x < xmin) {
2887 xmin = nodes[elements[imap].emap[6]].x;
2888 }
2889 if (nodes[elements[imap].emap[6]].x > xmax) {
2890 xmax = nodes[elements[imap].emap[6]].x;
2891 }
2892 if (nodes[elements[imap].emap[7]].x < xmin) {
2893 xmin = nodes[elements[imap].emap[7]].x;
2894 }
2895 if (nodes[elements[imap].emap[7]].x > xmax) {
2896 xmax = nodes[elements[imap].emap[7]].x;
2897 }
2898 ymin = nodes[elements[imap].emap[0]].y;
2899 ymax = nodes[elements[imap].emap[0]].y;
2900 if (nodes[elements[imap].emap[1]].y < ymin) {
2901 ymin = nodes[elements[imap].emap[1]].y;
2902 }
2903 if (nodes[elements[imap].emap[1]].y > ymax) {
2904 ymax = nodes[elements[imap].emap[1]].y;
2905 }
2906 if (nodes[elements[imap].emap[2]].y < ymin) {
2907 ymin = nodes[elements[imap].emap[2]].y;
2908 }
2909 if (nodes[elements[imap].emap[2]].y > ymax) {
2910 ymax = nodes[elements[imap].emap[2]].y;
2911 }
2912 if (nodes[elements[imap].emap[3]].y < ymin) {
2913 ymin = nodes[elements[imap].emap[3]].y;
2914 }
2915 if (nodes[elements[imap].emap[3]].y > ymax) {
2916 ymax = nodes[elements[imap].emap[3]].y;
2917 }
2918 if (nodes[elements[imap].emap[4]].y < ymin) {
2919 ymin = nodes[elements[imap].emap[4]].y;
2920 }
2921 if (nodes[elements[imap].emap[4]].y > ymax) {
2922 ymax = nodes[elements[imap].emap[4]].y;
2923 }
2924 if (nodes[elements[imap].emap[5]].y < ymin) {
2925 ymin = nodes[elements[imap].emap[5]].y;
2926 }
2927 if (nodes[elements[imap].emap[5]].y > ymax) {
2928 ymax = nodes[elements[imap].emap[5]].y;
2929 }
2930 if (nodes[elements[imap].emap[6]].y < ymin) {
2931 ymin = nodes[elements[imap].emap[6]].y;
2932 }
2933 if (nodes[elements[imap].emap[6]].y > ymax) {
2934 ymax = nodes[elements[imap].emap[6]].y;
2935 }
2936 if (nodes[elements[imap].emap[7]].y < ymin) {
2937 ymin = nodes[elements[imap].emap[7]].y;
2938 }
2939 if (nodes[elements[imap].emap[7]].y > ymax) {
2940 ymax = nodes[elements[imap].emap[7]].y;
2941 }
2942
2943 if (x >= xmin && x <= xmax && y >= ymin && y <= ymax) {
2944 std::cout << m_className << "::Coordinates5:\n";
2945 std::cout << " No convergence achieved "
2946 << "when refining internal isoparametric coordinates\n";
2947 std::cout << " in element " << imap << " at position (" << x << ", "
2948 << y << ").\n";
2949 t1 = t2 = 0;
2950 return ifail;
2951 }
2952 }
2953
2954 // Convergence reached.
2955 t1 = td1;
2956 t2 = td2;
2957 t3 = 0;
2958 t4 = 0;
2959 if (debug) {
2960 std::cout << m_className << "::Coordinates5:\n";
2961 std::cout << " Convergence reached at (t1, t2) = (" << t1 << ", " << t2
2962 << ").\n";
2963 }
2964
2965 // For debugging purposes, show position.
2966 if (debug) {
2967 double xr =
2968 nodes[elements[imap].emap[0]].x *
2969 (-(1 - t1) * (1 - t2) * (1 + t1 + t2)) / 4 +
2970 nodes[elements[imap].emap[1]].x *
2971 (-(1 + t1) * (1 - t2) * (1 - t1 + t2)) / 4 +
2972 nodes[elements[imap].emap[2]].x *
2973 (-(1 + t1) * (1 + t2) * (1 - t1 - t2)) / 4 +
2974 nodes[elements[imap].emap[3]].x *
2975 (-(1 - t1) * (1 + t2) * (1 + t1 - t2)) / 4 +
2976 nodes[elements[imap].emap[4]].x * (1 - t1) * (1 + t1) * (1 - t2) / 2 +
2977 nodes[elements[imap].emap[5]].x * (1 + t1) * (1 + t2) * (1 - t2) / 2 +
2978 nodes[elements[imap].emap[6]].x * (1 - t1) * (1 + t1) * (1 + t2) / 2 +
2979 nodes[elements[imap].emap[7]].x * (1 - t1) * (1 + t2) * (1 - t2) / 2;
2980 double yr =
2981 nodes[elements[imap].emap[0]].y *
2982 (-(1 - t1) * (1 - t2) * (1 + t1 + t2)) / 4 +
2983 nodes[elements[imap].emap[1]].y *
2984 (-(1 + t1) * (1 - t2) * (1 - t1 + t2)) / 4 +
2985 nodes[elements[imap].emap[2]].y *
2986 (-(1 + t1) * (1 + t2) * (1 - t1 - t2)) / 4 +
2987 nodes[elements[imap].emap[3]].y *
2988 (-(1 - t1) * (1 + t2) * (1 + t1 - t2)) / 4 +
2989 nodes[elements[imap].emap[4]].y * (1 - t1) * (1 + t1) * (1 - t2) / 2 +
2990 nodes[elements[imap].emap[5]].y * (1 + t1) * (1 + t2) * (1 - t2) / 2 +
2991 nodes[elements[imap].emap[6]].y * (1 - t1) * (1 + t1) * (1 + t2) / 2 +
2992 nodes[elements[imap].emap[7]].y * (1 - t1) * (1 + t2) * (1 - t2) / 2;
2993 std::cout << m_className << "::Coordinates5:\n";
2994 std::cout << " Position requested: (" << x << ", " << y << ")\n";
2995 std::cout << " Reconstructed: (" << xr << ", " << yr << ")\n";
2996 std::cout << " Difference: (" << x - xr << ", " << y - yr
2997 << ")\n";
2998 }
2999
3000 // Success
3001 ifail = 0;
3002 return ifail;
3003}
Garfield::ComponentFieldMap::Jacobian5
void Jacobian5(int i, double u, double v, double &det, double jac[4][4])
Definition: ComponentFieldMap.cc:696
Garfield::ComponentFieldMap::Coordinates4
int Coordinates4(double x, double y, double z, double &t1, double &t2, double &t3, double &t4, double jac[4][4], double &det, int imap)
Definition: ComponentFieldMap.cc:2460

Referenced by FindElement5().

◆ CoordinatesCube()

int Garfield::ComponentFieldMap::CoordinatesCube ( double  x,
double  y,
double  z,
double &  t1,
double &  t2,
double &  t3,
TMatrixD *&  jac,
std::vector< TMatrixD * > &  dN,
int  imap 
)
protected

Definition at line 3528 of file ComponentFieldMap.cc.

3530 {
3531
3532 /*
3533 global coordinates 7__ _ _ 6 t3 t2
3534 / /| ^ /|
3535 ^ z / / | | /
3536 | 4_______5 | | /
3537 | | | | | /
3538 | | 3 | 2 |/ t1
3539 -------> | | / ------->
3540 / y | |/ local coordinates
3541 / 0--------1
3542 /
3543 v x
3544 */
3545
3546 // Failure flag
3547 int ifail = 1;
3548
3549 // Compute hexahedral coordinates (t1->[-1,1],t2->[-1,1],t3->[-1,1]) and
3550 // t1 (zeta) is in y-direction
3551 // t2 (eta) is in opposite x-direction
3552 // t3 (mu) is in z-direction
3553 // Nodes are set in that way, that node [0] has always lowest x,y,z!
3554 t2 =
3555 (2. * (x - nodes[elements[imap].emap[3]].x) /
3556 (nodes[elements[imap].emap[0]].x - nodes[elements[imap].emap[3]].x) -
3557 1) *
3558 -1.;
3559 t1 = 2. * (y - nodes[elements[imap].emap[3]].y) /
3560 (nodes[elements[imap].emap[2]].y - nodes[elements[imap].emap[3]].y) -
3561 1;
3562 t3 = 2. * (z - nodes[elements[imap].emap[3]].z) /
3563 (nodes[elements[imap].emap[7]].z - nodes[elements[imap].emap[3]].z) -
3564 1;
3565 // Re-compute the (x,y,z) position for this coordinate.
3566 if (debug) {
3567 double n[8];
3568 n[0] = 1. / 8 * (1 - t1) * (1 - t2) * (1 - t3);
3569 n[1] = 1. / 8 * (1 + t1) * (1 - t2) * (1 - t3);
3570 n[2] = 1. / 8 * (1 + t1) * (1 + t2) * (1 - t3);
3571 n[3] = 1. / 8 * (1 - t1) * (1 + t2) * (1 - t3);
3572 n[4] = 1. / 8 * (1 - t1) * (1 - t2) * (1 + t3);
3573 n[5] = 1. / 8 * (1 + t1) * (1 - t2) * (1 + t3);
3574 n[6] = 1. / 8 * (1 + t1) * (1 + t2) * (1 + t3);
3575 n[7] = 1. / 8 * (1 - t1) * (1 + t2) * (1 + t3);
3576
3577 double xr = 0;
3578 double yr = 0;
3579 double zr = 0;
3580
3581 for (int i = 0; i < 8; i++) {
3582 xr += nodes[elements[imap].emap[i]].x * n[i];
3583 yr += nodes[elements[imap].emap[i]].y * n[i];
3584 zr += nodes[elements[imap].emap[i]].z * n[i];
3585 }
3586 double sr = n[0] + n[1] + n[2] + n[3] + n[4] + n[5] + n[6] + n[7];
3587 std::cout << m_className << "::CoordinatesCube:\n";
3588 std::cout << " Position requested: (" << x << "," << y << "," << z
3589 << ")\n";
3590 std::cout << " Position reconstructed: (" << xr << "," << yr << "," << zr
3591 << ")\n";
3592 std::cout << " Difference: (" << (x - xr) << "," << (y - yr)
3593 << "," << (z - zr) << ")\n";
3594 std::cout << " Hexahedral coordinates (t, u, v) = (" << t1 << "," << t2
3595 << "," << t3 << ")\n";
3596 std::cout << " Checksum - 1: " << (sr - 1) << "\n";
3597 }
3598 if (jac != 0) JacobianCube(imap, t1, t2, t3, jac, dN);
3599 // This should always work.
3600 ifail = 0;
3601 return ifail;
3602}
Garfield::ComponentFieldMap::JacobianCube
void JacobianCube(int i, double t1, double t2, double t3, TMatrixD *&jac, std::vector< TMatrixD * > &dN)
Definition: ComponentFieldMap.cc:2128

Referenced by FindElementCube().

◆ DisableCheckMapIndices()

void Garfield::ComponentFieldMap::DisableCheckMapIndices ( )
inline

Definition at line 76 of file ComponentFieldMap.hh.

76{ checkMultipleElement = false; }

◆ DisableDeleteBackgroundElements()

void Garfield::ComponentFieldMap::DisableDeleteBackgroundElements ( )
inline

Definition at line 78 of file ComponentFieldMap.hh.

78{ deleteBackground = false; }

◆ DriftMedium()

void Garfield::ComponentFieldMap::DriftMedium ( int  imat)

Definition at line 71 of file ComponentFieldMap.cc.

71 {
72
73 // Do not proceed if not properly initialised.
74 if (!ready) {
75 std::cerr << m_className << "::DriftMedium:\n";
76 std::cerr << " Field map not yet initialised.\n";
77 std::cerr << " Drift medium cannot be selected.\n";
78 return;
79 }
80
81 // Check value
82 if (imat < 0 || imat >= nMaterials) {
83 std::cerr << m_className << "::DriftMedium:\n";
84 std::cerr << " Material index " << imat << " is out of range.\n";
85 return;
86 }
87
88 // Make drift medium
89 materials[imat].driftmedium = true;
90}

◆ ElectricField() [1/2]

virtual void Garfield::ComponentFieldMap::ElectricField ( const double  x,
const double  y,
const double  z,
double &  ex,
double &  ey,
double &  ez,
double &  v,
Medium *&  m,
int &  status 
)
pure virtual

Implements Garfield::ComponentBase.

Implemented in Garfield::ComponentAnsys121, Garfield::ComponentAnsys123, Garfield::ComponentCST, and Garfield::ComponentElmer.

◆ ElectricField() [2/2]

virtual void Garfield::ComponentFieldMap::ElectricField ( const double  x,
const double  y,
const double  z,
double &  ex,
double &  ey,
double &  ez,
Medium *&  m,
int &  status 
)
pure virtual

Implements Garfield::ComponentBase.

Implemented in Garfield::ComponentAnsys121, Garfield::ComponentAnsys123, Garfield::ComponentCST, and Garfield::ComponentElmer.

◆ EnableCheckMapIndices()

void Garfield::ComponentFieldMap::EnableCheckMapIndices ( )
inline

Definition at line 72 of file ComponentFieldMap.hh.

72 {
73 checkMultipleElement = true;
74 lastElement = -1;
75 }

◆ EnableDeleteBackgroundElements()

void Garfield::ComponentFieldMap::EnableDeleteBackgroundElements ( )
inline

Definition at line 77 of file ComponentFieldMap.hh.

77{ deleteBackground = true; }

◆ FindElement13()

int Garfield::ComponentFieldMap::FindElement13 ( const double  x,
const double  y,
const double  z,
double &  t1,
double &  t2,
double &  t3,
double &  t4,
double  jac[4][4],
double &  det 
)
protected

Definition at line 356 of file ComponentFieldMap.cc.

359 {
360 // Check if bounding boxes of elements have been computed
361 if(!cacheElemBoundingBoxes) {
362 std::cerr << m_className << "::FindElement13:\n";
363 std::cerr << " Caching the bounding box calculations of all elements.\n";
364
365 CalculateElementBoundingBoxes();
366 cacheElemBoundingBoxes = true;
367 }
368
369 // Backup
370 double jacbak[4][4];
371 double detbak = 1.;
372 double t1bak = 0., t2bak = 0., t3bak = 0., t4bak = 0.;
373 int imapbak = -1;
374
375 // Initial values.
376 t1 = t2 = t3 = t4 = 0.;
377
378 // Check previously used element
379 int rc;
380 if (lastElement > -1 && !checkMultipleElement) {
381 rc = Coordinates13(x, y, z, t1, t2, t3, t4, jac, det, lastElement);
382 if (rc == 0 && t1 >= 0 && t1 <= +1 && t2 >= 0 && t2 <= +1 && t3 >= 0 &&
383 t3 <= +1 && t4 >= 0 && t4 <= +1)
384 return lastElement;
385 }
386
387 // Verify the count of volumes that contain the point.
388 int nfound = 0;
389 int imap = -1;
390
391 // Tolerance
392 const double f = 0.2;
393
394 // Scan all elements
395 for (int i = 0; i < nElements; i++) {
396 element& e = elements[i];
397 if (x < e.xmin - f * (e.xmax - e.xmin) || x > e.xmax + f * (e.xmax - e.xmin) ||
398 y < e.ymin - f * (e.ymax - e.ymin) || y > e.ymax + f * (e.ymax - e.ymin) ||
399 z < e.zmin - f * (e.zmax - e.zmin) || z > e.zmax + f * (e.zmax - e.zmin))
400 continue;
401
402 rc = Coordinates13(x, y, z, t1, t2, t3, t4, jac, det, i);
403
404 if (rc == 0 && t1 >= 0 && t1 <= +1 && t2 >= 0 && t2 <= +1 && t3 >= 0 &&
405 t3 <= +1 && t4 >= 0 && t4 <= +1) {
406 ++nfound;
407 imap = i;
408 lastElement = i;
409 if (debug) {
410 std::cout << m_className << "::FindElement13:\n";
411 std::cout << " Found matching element " << i << ".\n";
412 }
413 if (!checkMultipleElement) return i;
414 for (int j = 0; j < 4; ++j) {
415 for (int k = 0; k < 4; ++k) jacbak[j][k] = jac[j][k];
416 }
417 detbak = det;
418 t1bak = t1;
419 t2bak = t2;
420 t3bak = t3;
421 t4bak = t4;
422 imapbak = imap;
423 if (debug) {
424 std::cout << m_className << "::FindElement13:\n";
425 std::cout << " Global = (" << x << ", " << y << ")\n";
426 std::cout << " Local = (" << t1 << ", " << t2 << ", " << t3 << ", "
427 << t4 << ")\n";
428 std::cout << " Element = " << imap << "\n";
429 std::cout << " Node x "
430 "y z V\n";
431 for (int ii = 0; ii < 10; ++ii) {
432 printf(" %-5d %12g %12g %12g %12g\n",
433 elements[imap].emap[ii], nodes[elements[imap].emap[ii]].x,
434 nodes[elements[imap].emap[ii]].y,
435 nodes[elements[imap].emap[ii]].z,
436 nodes[elements[imap].emap[ii]].v);
437 }
438 }
439 }
440 }
441
442 // In checking mode, verify the tetrahedron/triangle count.
443 if (checkMultipleElement) {
444 if (nfound < 1) {
445 if (debug) {
446 std::cout << m_className << "::FindElement13:\n";
447 std::cout << " No element matching point (" << x << ", " << y << ", "
448 << z << ") found.\n";
449 }
450 lastElement = -1;
451 return -1;
452 }
453 if (nfound > 1) {
454 std::cerr << m_className << "::FindElement13:\n";
455 std::cerr << " Found << " << nfound << " elements matching point ("
456 << x << ", " << y << ", " << z << ").\n";
457 }
458 if (nfound > 0) {
459 for (int j = 0; j < 4; ++j) {
460 for (int k = 0; k < 4; ++k) jac[j][k] = jacbak[j][k];
461 }
462 det = detbak;
463 t1 = t1bak;
464 t2 = t2bak;
465 t3 = t3bak;
466 t4 = t4bak;
467 imap = imapbak;
468 lastElement = imap;
469 return imap;
470 }
471 }
472
473 if (debug) {
474 std::cout << m_className << "::FindElement13:\n";
475 std::cout << " No element matching point (" << x << ", " << y << ", "
476 << z << ") found.\n";
477 }
478 return -1;
479}
Garfield::ComponentFieldMap::Coordinates13
int Coordinates13(double x, double y, double z, double &t1, double &t2, double &t3, double &t4, double jac[4][4], double &det, int imap)
Definition: ComponentFieldMap.cc:3257

Referenced by Garfield::ComponentAnsys123::ElectricField(), Garfield::ComponentElmer::ElectricField(), Garfield::ComponentAnsys123::GetMedium(), Garfield::ComponentElmer::GetMedium(), Garfield::ComponentAnsys123::WeightingField(), Garfield::ComponentElmer::WeightingField(), Garfield::ComponentAnsys123::WeightingPotential(), and Garfield::ComponentElmer::WeightingPotential().

◆ FindElement5()

int Garfield::ComponentFieldMap::FindElement5 ( const double  x,
const double  y,
const double  z,
double &  t1,
double &  t2,
double &  t3,
double &  t4,
double  jac[4][4],
double &  det 
)
protected

Definition at line 184 of file ComponentFieldMap.cc.

187 {
188
189 // Check if bounding boxes of elements have been computed
190 if(!cacheElemBoundingBoxes) {
191 std::cerr << m_className << "::FindElement5:\n";
192 std::cerr << " Caching the bounding box calculations of all elements.\n";
193
194 CalculateElementBoundingBoxes();
195 cacheElemBoundingBoxes = true;
196 }
197
198 // Backup
199 double jacbak[4][4], detbak = 1.;
200 double t1bak = 0., t2bak = 0., t3bak = 0., t4bak = 0.;
201 int imapbak = -1;
202
203 // Initial values.
204 t1 = t2 = t3 = t4 = 0;
205
206 // Check previously used element
207 int rc;
208 if (lastElement > -1 && !checkMultipleElement) {
209 if (elements[lastElement].degenerate) {
210 rc = Coordinates3(x, y, z, t1, t2, t3, t4, jac, det, lastElement);
211 if (rc == 0 && t1 >= 0 && t1 <= +1 && t2 >= 0 && t2 <= +1 && t3 >= 0 &&
212 t3 <= +1)
213 return lastElement;
214 } else {
215 rc = Coordinates5(x, y, z, t1, t2, t3, t4, jac, det, lastElement);
216 if (rc == 0 && t1 >= -1 && t1 <= +1 && t2 >= -1 && t2 <= +1)
217 return lastElement;
218 }
219 }
220
221 // Verify the count of volumes that contain the point.
222 int nfound = 0;
223 int imap = -1;
224
225 // Tolerance
226 const double f = 0.2;
227
228 // Scan all elements
229 for (int i = 0; i < nElements; ++i) {
230 element& e = elements[i];
231 if (x < e.xmin - f * (e.xmax - e.xmin) || x > e.xmax + f * (e.xmax - e.xmin) ||
232 y < e.ymin - f * (e.ymax - e.ymin) || y > e.ymax + f * (e.ymax - e.ymin) ||
233 z < e.zmin - f * (e.zmax - e.zmin) || z > e.zmax + f * (e.zmax - e.zmin))
234 continue;
235
236 if (elements[i].degenerate) {
237 // Degenerate element
238 rc = Coordinates3(x, y, z, t1, t2, t3, t4, jac, det, i);
239 if (rc == 0 && t1 >= 0 && t1 <= +1 && t2 >= 0 && t2 <= +1 && t3 >= 0 &&
240 t3 <= +1) {
241 ++nfound;
242 imap = i;
243 lastElement = i;
244 if (debug) {
245 std::cout << m_className << "::FindElement5:\n";
246 std::cout << " Found matching degenerate element " << i << ".\n";
247 }
248 if (!checkMultipleElement) return i;
249 for (int j = 0; j < 4; ++j) {
250 for (int k = 0; k < 4; ++k) jacbak[j][k] = jac[j][k];
251 }
252 detbak = det;
253 t1bak = t1;
254 t2bak = t2;
255 t3bak = t3;
256 t4bak = t4;
257 imapbak = imap;
258 if (debug) {
259 std::cout << m_className << "::FindElement5:\n";
260 std::cout << " Global = (" << x << ", " << y << ")\n";
261 std::cout << " Local = (" << t1 << ", " << t2 << ", " << t3 << ", "
262 << t4 << ")\n";
263 std::cout << " Element = " << imap
264 << " (degenerate: " << elements[imap].degenerate << ")\n";
265 std::cout << " Node x "
266 " y V\n";
267 for (int ii = 0; ii < 6; ++ii) {
268 printf(" %-5d %12g %12g %12g\n",
269 elements[imap].emap[ii], nodes[elements[imap].emap[ii]].x,
270 nodes[elements[imap].emap[ii]].y,
271 nodes[elements[imap].emap[ii]].v);
272 }
273 }
274 }
275 } else {
276 // Non-degenerate element
277 rc = Coordinates5(x, y, z, t1, t2, t3, t4, jac, det, i);
278 if (rc == 0 && t1 >= -1 && t1 <= +1 && t2 >= -1 && t2 <= +1) {
279 ++nfound;
280 imap = i;
281 lastElement = i;
282 if (debug) {
283 std::cout << m_className << "::FindElement5:\n";
284 std::cout << " Found matching non-degenerate element " << i
285 << ".\n";
286 }
287 if (!checkMultipleElement) return i;
288 for (int j = 0; j < 4; ++j) {
289 for (int k = 0; k < 4; ++k) jacbak[j][k] = jac[j][k];
290 }
291 detbak = det;
292 t1bak = t1;
293 t2bak = t2;
294 t3bak = t3;
295 t4bak = t4;
296 imapbak = imap;
297 if (debug) {
298 std::cout << m_className << "::FindElement5:\n";
299 std::cout << " Global = (" << x << ", " << y << ")\n";
300 std::cout << " Local = (" << t1 << ", " << t2 << ", " << t3 << ", "
301 << t4 << ")\n";
302 std::cout << " Element = " << imap
303 << " (degenerate: " << elements[imap].degenerate << ")\n";
304 std::cout << " Node x "
305 " y V\n";
306 for (int ii = 0; ii < 8; ++ii) {
307 printf(" %-5d %12g %12g %12g\n",
308 elements[imap].emap[ii], nodes[elements[imap].emap[ii]].x,
309 nodes[elements[imap].emap[ii]].y,
310 nodes[elements[imap].emap[ii]].v);
311 }
312 }
313 }
314 }
315 }
316
317 // In checking mode, verify the tetrahedron/triangle count.
318 if (checkMultipleElement) {
319 if (nfound < 1) {
320 if (debug) {
321 std::cout << m_className << "::FindElement5:\n";
322 std::cout << " No element matching point (" << x << ", " << y
323 << ") found.\n";
324 }
325 lastElement = -1;
326 return -1;
327 }
328 if (nfound > 1) {
329 std::cout << m_className << "::FindElement5:\n";
330 std::cout << " Found " << nfound << " elements matching point (" << x
331 << ", " << y << ").\n";
332 }
333 if (nfound > 0) {
334 for (int j = 0; j < 4; ++j) {
335 for (int k = 0; k < 4; ++k) jac[j][k] = jacbak[j][k];
336 }
337 det = detbak;
338 t1 = t1bak;
339 t2 = t2bak;
340 t3 = t3bak;
341 t4 = t4bak;
342 imap = imapbak;
343 lastElement = imap;
344 return imap;
345 }
346 }
347
348 if (debug) {
349 std::cout << m_className << "::FindElement5:\n";
350 std::cout << " No element matching point (" << x << ", " << y
351 << ") found.\n";
352 }
353 return -1;
354}
Garfield::ComponentFieldMap::Coordinates5
int Coordinates5(double x, double y, double z, double &t1, double &t2, double &t3, double &t4, double jac[4][4], double &det, int imap)
Definition: ComponentFieldMap.cc:2721
Garfield::ComponentFieldMap::Coordinates3
int Coordinates3(double x, double y, double z, double &t1, double &t2, double &t3, double &t4, double jac[4][4], double &det, int imap)
Definition: ComponentFieldMap.cc:2226

Referenced by Garfield::ComponentAnsys121::ElectricField(), Garfield::ComponentAnsys121::GetMedium(), Garfield::ComponentAnsys121::WeightingField(), and Garfield::ComponentAnsys121::WeightingPotential().

◆ FindElementCube()

int Garfield::ComponentFieldMap::FindElementCube ( const double  x,
const double  y,
const double  z,
double &  t1,
double &  t2,
double &  t3,
TMatrixD *&  jac,
std::vector< TMatrixD * > &  dN 
)
protected

Definition at line 481 of file ComponentFieldMap.cc.

484 {
485
486 int imap = -1;
487
488 if (x >= nodes[elements[lastElement].emap[3]].x &&
489 y >= nodes[elements[lastElement].emap[3]].y &&
490 z >= nodes[elements[lastElement].emap[3]].z &&
491 x < nodes[elements[lastElement].emap[0]].x &&
492 y < nodes[elements[lastElement].emap[2]].y &&
493 z < nodes[elements[lastElement].emap[7]].z) {
494 imap = lastElement;
495 }
496
497 // Default element loop
498 if (imap == -1) {
499 for (int i = 0; i < nElements; ++i) {
500 if (x >= nodes[elements[i].emap[3]].x &&
501 y >= nodes[elements[i].emap[3]].y &&
502 z >= nodes[elements[i].emap[3]].z &&
503 x < nodes[elements[i].emap[0]].x &&
504 y < nodes[elements[i].emap[2]].y &&
505 z < nodes[elements[i].emap[7]].z) {
506 imap = i;
507 break;
508 }
509 }
510 }
511
512 if (imap < 0) {
513 if (debug) {
514 std::cout << m_className << "::FindElementCube:\n";
515 std::cout << " Point (" << x << "," << y << "," << z
516 << ") not in the mesh, it is background or PEC.\n";
517 std::cout << " First node (" << nodes[elements[0].emap[3]].x << ","
518 << nodes[elements[0].emap[3]].y << ","
519 << nodes[elements[0].emap[3]].z << ") in the mesh.\n";
520 std::cout << " dx= "
521 << (nodes[elements[0].emap[0]].x - nodes[elements[0].emap[3]].x)
522 << ", dy= "
523 << (nodes[elements[0].emap[2]].y - nodes[elements[0].emap[3]].y)
524 << ", dz= "
525 << (nodes[elements[0].emap[7]].z - nodes[elements[0].emap[3]].z)
526 << "\n";
527 std::cout << " Last node (" << nodes[elements[nElements - 1].emap[5]].x
528 << "," << nodes[elements[nElements - 1].emap[5]].y << ","
529 << nodes[elements[nElements - 1].emap[5]].z
530 << ") in the mesh.\n";
531 std::cout << " dx= " << (nodes[elements[nElements - 1].emap[0]].x -
532 nodes[elements[nElements - 1].emap[3]].x)
533 << ", dy= " << (nodes[elements[nElements - 1].emap[2]].y -
534 nodes[elements[nElements - 1].emap[3]].y)
535 << ", dz= " << (nodes[elements[nElements - 1].emap[7]].z -
536 nodes[elements[nElements - 1].emap[3]].z)
537 << "\n";
538 }
539 return -1;
540 }
541 CoordinatesCube(x, y, z, t1, t2, t3, jac, dN, imap);
542 if (debug) {
543 std::cout << m_className << "::FindElementCube:\n";
544 std::cout << "Global: (" << x << "," << y << "," << z << ") in element "
545 << imap << " (degenerate: " << elements[imap].degenerate << ")\n";
546 std::cout << " Node xyzV\n";
547 for (int i = 0; i < 8; i++) {
548 std::cout << " " << elements[imap].emap[i] << " "
549 << nodes[elements[imap].emap[i]].x << " "
550 << nodes[elements[imap].emap[i]].y << " "
551 << nodes[elements[imap].emap[i]].z << " "
552 << nodes[elements[imap].emap[i]].v << "\n";
553 }
554 }
555 return imap;
556}
Garfield::ComponentFieldMap::CoordinatesCube
int CoordinatesCube(double x, double y, double z, double &t1, double &t2, double &t3, TMatrixD *&jac, std::vector< TMatrixD * > &dN, int imap)
Definition: ComponentFieldMap.cc:3528

◆ GetAspectRatio()

virtual void Garfield::ComponentFieldMap::GetAspectRatio ( const int  i,
double &  dmin,
double &  dmax 
)
protectedpure virtual

Implemented in Garfield::ComponentAnsys121, Garfield::ComponentAnsys123, Garfield::ComponentCST, and Garfield::ComponentElmer.

Referenced by GetElement().

◆ GetBoundingBox()

bool Garfield::ComponentFieldMap::GetBoundingBox ( double &  xmin,
double &  ymin,
double &  zmin,
double &  xmax,
double &  ymax,
double &  zmax 
)
virtual

Reimplemented from Garfield::ComponentBase.

Definition at line 4021 of file ComponentFieldMap.cc.

4023 {
4024
4025 if (!ready) return false;
4026
4027 xmin = xMinBoundingBox;
4028 xmax = xMaxBoundingBox;
4029 ymin = yMinBoundingBox;
4030 ymax = yMaxBoundingBox;
4031 zmin = zMinBoundingBox;
4032 zmax = zMaxBoundingBox;
4033 return true;
4034}

◆ GetConductivity()

double Garfield::ComponentFieldMap::GetConductivity ( const int  imat)

Definition at line 124 of file ComponentFieldMap.cc.

124 {
125
126 if (imat < 0 || imat >= nMaterials) {
127 std::cerr << m_className << "::GetConductivity:\n";
128 std::cerr << " Material index " << imat << " is out of range.\n";
129 return -1.;
130 }
131
132 return materials[imat].ohm;
133}

◆ GetElement()

bool Garfield::ComponentFieldMap::GetElement ( const int  i,
double &  vol,
double &  dmin,
double &  dmax 
)

Definition at line 170 of file ComponentFieldMap.cc.

171 {
172
173 if (i < 0 || i >= nElements) {
174 std::cerr << m_className << "::GetElement:\n";
175 std::cerr << " Element index (" << i << ") out of range.\n";
176 return false;
177 }
178
179 vol = GetElementVolume(i);
180 GetAspectRatio(i, dmin, dmax);
181 return true;
182}
Garfield::ComponentFieldMap::GetElementVolume
virtual double GetElementVolume(const int i)=0
Garfield::ComponentFieldMap::GetAspectRatio
virtual void GetAspectRatio(const int i, double &dmin, double &dmax)=0

◆ GetElementVolume()

virtual double Garfield::ComponentFieldMap::GetElementVolume ( const int  i)
protectedpure virtual

Implemented in Garfield::ComponentAnsys121, Garfield::ComponentAnsys123, Garfield::ComponentCST, and Garfield::ComponentElmer.

Referenced by GetElement().

◆ GetMedium() [1/2]

Medium * Garfield::ComponentFieldMap::GetMedium ( const double &  x,
const double &  y,
const double &  z 
)
pure virtual

Reimplemented from Garfield::ComponentBase.

Implemented in Garfield::ComponentAnsys121, Garfield::ComponentAnsys123, Garfield::ComponentCST, and Garfield::ComponentElmer.

◆ GetMedium() [2/2]

Medium * Garfield::ComponentFieldMap::GetMedium ( const unsigned int &  i) const

Definition at line 159 of file ComponentFieldMap.cc.

159 {
160
161 if (imat >= (unsigned int)nMaterials) {
162 std::cerr << m_className << "::GetMedium:\n";
163 std::cerr << " Material index " << imat << " is out of range.\n";
164 return NULL;
165 }
166
167 return materials[imat].medium;
168}

◆ GetNumberOfElements()

int Garfield::ComponentFieldMap::GetNumberOfElements ( ) const
inline

Definition at line 53 of file ComponentFieldMap.hh.

53{ return nElements; }

◆ GetNumberOfMaterials()

int Garfield::ComponentFieldMap::GetNumberOfMaterials ( )
inline

Definition at line 41 of file ComponentFieldMap.hh.

41{ return nMaterials; }

◆ GetNumberOfMedia()

int Garfield::ComponentFieldMap::GetNumberOfMedia ( )
inline

Definition at line 51 of file ComponentFieldMap.hh.

51{ return nMaterials; }

◆ GetPermittivity()

double Garfield::ComponentFieldMap::GetPermittivity ( const int  imat)

Definition at line 113 of file ComponentFieldMap.cc.

113 {
114
115 if (imat < 0 || imat >= nMaterials) {
116 std::cerr << m_className << "::GetPermittivity:\n";
117 std::cerr << " Material index " << imat << " is out of range.\n";
118 return -1.;
119 }
120
121 return materials[imat].eps;
122}

◆ GetVoltageRange()

bool Garfield::ComponentFieldMap::GetVoltageRange ( double &  vmin,
double &  vmax 
)
inlinevirtual

Implements Garfield::ComponentBase.

Definition at line 27 of file ComponentFieldMap.hh.

27 {
28 vmin = mapvmin;
29 vmax = mapvmax;
30 return true;
31 }

Referenced by Garfield::ViewFEMesh::Plot().

◆ IsInBoundingBox()

bool Garfield::ComponentFieldMap::IsInBoundingBox ( const double  x,
const double  y,
const double  z 
)
virtual

Reimplemented in Garfield::ComponentAnsys121, Garfield::ComponentAnsys123, Garfield::ComponentCST, and Garfield::ComponentElmer.

Definition at line 4011 of file ComponentFieldMap.cc.

4012 {
4013
4014 if (x >= xMinBoundingBox && x <= xMaxBoundingBox && y >= yMinBoundingBox &&
4015 y <= yMaxBoundingBox && z >= zMinBoundingBox && z <= zMaxBoundingBox) {
4016 return true;
4017 }
4018 return false;
4019}

◆ Jacobian13()

void Garfield::ComponentFieldMap::Jacobian13 ( int  i,
double  t,
double  u,
double  v,
double  w,
double &  det,
double  jac[4][4] 
)
protected

Definition at line 1233 of file ComponentFieldMap.cc.

1234 {
1235
1236 // Initial values
1237 det = 0;
1238 for (int j = 0; j < 4; ++j) {
1239 for (int k = 0; k < 4; ++k) jac[j][k] = 0;
1240 }
1241
1242 // Be sure that the element is within range
1243 if (i < 0 || i >= nElements) {
1244 std::cerr << m_className << "::Jacobian13:\n";
1245 std::cerr << " Element " << i << " out of range.\n";
1246 return;
1247 }
1248
1249 // Determinant of the quadrilateral serendipity Jacobian
1250 det =
1251 -(((-4 * v * nodes[elements[i].emap[9]].x - nodes[elements[i].emap[1]].x +
1252 4 * u * nodes[elements[i].emap[1]].x + nodes[elements[i].emap[3]].x -
1253 4 * w * nodes[elements[i].emap[3]].x +
1254 4 * t * nodes[elements[i].emap[4]].x -
1255 4 * t * nodes[elements[i].emap[6]].x +
1256 4 * v * nodes[elements[i].emap[7]].x -
1257 4 * u * nodes[elements[i].emap[8]].x +
1258 4 * w * nodes[elements[i].emap[8]].x) *
1259 (4 * w * nodes[elements[i].emap[9]].y -
1260 nodes[elements[i].emap[2]].y +
1261 4 * v * nodes[elements[i].emap[2]].y +
1262 4 * t * nodes[elements[i].emap[5]].y +
1263 4 * u * nodes[elements[i].emap[7]].y) +
1264 (nodes[elements[i].emap[1]].x - 4 * u * nodes[elements[i].emap[1]].x -
1265 nodes[elements[i].emap[2]].x + 4 * v * nodes[elements[i].emap[2]].x -
1266 4 * t * nodes[elements[i].emap[4]].x +
1267 4 * t * nodes[elements[i].emap[5]].x +
1268 4 * u * nodes[elements[i].emap[7]].x -
1269 4 * v * nodes[elements[i].emap[7]].x +
1270 4 * w *
1271 (nodes[elements[i].emap[9]].x - nodes[elements[i].emap[8]].x)) *
1272 (4 * v * nodes[elements[i].emap[9]].y -
1273 nodes[elements[i].emap[3]].y +
1274 4 * w * nodes[elements[i].emap[3]].y +
1275 4 * t * nodes[elements[i].emap[6]].y +
1276 4 * u * nodes[elements[i].emap[8]].y) +
1277 (-4 * w * nodes[elements[i].emap[9]].x +
1278 4 * v *
1279 (nodes[elements[i].emap[9]].x - nodes[elements[i].emap[2]].x) +
1280 nodes[elements[i].emap[2]].x - nodes[elements[i].emap[3]].x +
1281 4 * w * nodes[elements[i].emap[3]].x -
1282 4 * t * nodes[elements[i].emap[5]].x +
1283 4 * t * nodes[elements[i].emap[6]].x -
1284 4 * u * nodes[elements[i].emap[7]].x +
1285 4 * u * nodes[elements[i].emap[8]].x) *
1286 ((-1 + 4 * u) * nodes[elements[i].emap[1]].y +
1287 4 * (t * nodes[elements[i].emap[4]].y +
1288 v * nodes[elements[i].emap[7]].y +
1289 w * nodes[elements[i].emap[8]].y))) *
1290 ((-1 + 4 * t) * nodes[elements[i].emap[0]].z +
1291 4 * (u * nodes[elements[i].emap[4]].z +
1292 v * nodes[elements[i].emap[5]].z +
1293 w * nodes[elements[i].emap[6]].z))) -
1294 ((nodes[elements[i].emap[1]].x - 4 * u * nodes[elements[i].emap[1]].x -
1295 nodes[elements[i].emap[3]].x + 4 * w * nodes[elements[i].emap[3]].x -
1296 4 * t * nodes[elements[i].emap[4]].x +
1297 4 * t * nodes[elements[i].emap[6]].x +
1298 4 * v * (nodes[elements[i].emap[9]].x - nodes[elements[i].emap[7]].x) +
1299 4 * u * nodes[elements[i].emap[8]].x -
1300 4 * w * nodes[elements[i].emap[8]].x) *
1301 ((-1 + 4 * t) * nodes[elements[i].emap[0]].y +
1302 4 * (u * nodes[elements[i].emap[4]].y +
1303 v * nodes[elements[i].emap[5]].y +
1304 w * nodes[elements[i].emap[6]].y)) -
1305 ((-1 + 4 * t) * nodes[elements[i].emap[0]].x +
1306 nodes[elements[i].emap[1]].x - 4 * u * nodes[elements[i].emap[1]].x +
1307 4 * (-(t * nodes[elements[i].emap[4]].x) +
1308 u * nodes[elements[i].emap[4]].x +
1309 v * nodes[elements[i].emap[5]].x +
1310 w * nodes[elements[i].emap[6]].x -
1311 v * nodes[elements[i].emap[7]].x -
1312 w * nodes[elements[i].emap[8]].x)) *
1313 (4 * v * nodes[elements[i].emap[9]].y -
1314 nodes[elements[i].emap[3]].y +
1315 4 * w * nodes[elements[i].emap[3]].y +
1316 4 * t * nodes[elements[i].emap[6]].y +
1317 4 * u * nodes[elements[i].emap[8]].y) +
1318 ((-1 + 4 * t) * nodes[elements[i].emap[0]].x -
1319 4 * v * nodes[elements[i].emap[9]].x + nodes[elements[i].emap[3]].x -
1320 4 * w * nodes[elements[i].emap[3]].x +
1321 4 * u * nodes[elements[i].emap[4]].x +
1322 4 * v * nodes[elements[i].emap[5]].x -
1323 4 * t * nodes[elements[i].emap[6]].x +
1324 4 * w * nodes[elements[i].emap[6]].x -
1325 4 * u * nodes[elements[i].emap[8]].x) *
1326 ((-1 + 4 * u) * nodes[elements[i].emap[1]].y +
1327 4 * (t * nodes[elements[i].emap[4]].y +
1328 v * nodes[elements[i].emap[7]].y +
1329 w * nodes[elements[i].emap[8]].y))) *
1330 (4 * w * nodes[elements[i].emap[9]].z - nodes[elements[i].emap[2]].z +
1331 4 * v * nodes[elements[i].emap[2]].z +
1332 4 * t * nodes[elements[i].emap[5]].z +
1333 4 * u * nodes[elements[i].emap[7]].z) +
1334 ((nodes[elements[i].emap[1]].x - 4 * u * nodes[elements[i].emap[1]].x -
1335 nodes[elements[i].emap[2]].x + 4 * v * nodes[elements[i].emap[2]].x -
1336 4 * t * nodes[elements[i].emap[4]].x +
1337 4 * t * nodes[elements[i].emap[5]].x +
1338 4 * u * nodes[elements[i].emap[7]].x -
1339 4 * v * nodes[elements[i].emap[7]].x +
1340 4 * w * (nodes[elements[i].emap[9]].x - nodes[elements[i].emap[8]].x)) *
1341 ((-1 + 4 * t) * nodes[elements[i].emap[0]].y +
1342 4 * (u * nodes[elements[i].emap[4]].y +
1343 v * nodes[elements[i].emap[5]].y +
1344 w * nodes[elements[i].emap[6]].y)) -
1345 ((-1 + 4 * t) * nodes[elements[i].emap[0]].x +
1346 nodes[elements[i].emap[1]].x - 4 * u * nodes[elements[i].emap[1]].x +
1347 4 * (-(t * nodes[elements[i].emap[4]].x) +
1348 u * nodes[elements[i].emap[4]].x +
1349 v * nodes[elements[i].emap[5]].x +
1350 w * nodes[elements[i].emap[6]].x -
1351 v * nodes[elements[i].emap[7]].x -
1352 w * nodes[elements[i].emap[8]].x)) *
1353 (4 * w * nodes[elements[i].emap[9]].y -
1354 nodes[elements[i].emap[2]].y +
1355 4 * v * nodes[elements[i].emap[2]].y +
1356 4 * t * nodes[elements[i].emap[5]].y +
1357 4 * u * nodes[elements[i].emap[7]].y) +
1358 ((-1 + 4 * t) * nodes[elements[i].emap[0]].x -
1359 4 * w * nodes[elements[i].emap[9]].x + nodes[elements[i].emap[2]].x -
1360 4 * v * nodes[elements[i].emap[2]].x +
1361 4 * u * nodes[elements[i].emap[4]].x -
1362 4 * t * nodes[elements[i].emap[5]].x +
1363 4 * v * nodes[elements[i].emap[5]].x +
1364 4 * w * nodes[elements[i].emap[6]].x -
1365 4 * u * nodes[elements[i].emap[7]].x) *
1366 ((-1 + 4 * u) * nodes[elements[i].emap[1]].y +
1367 4 * (t * nodes[elements[i].emap[4]].y +
1368 v * nodes[elements[i].emap[7]].y +
1369 w * nodes[elements[i].emap[8]].y))) *
1370 (4 * v * nodes[elements[i].emap[9]].z - nodes[elements[i].emap[3]].z +
1371 4 * w * nodes[elements[i].emap[3]].z +
1372 4 * t * nodes[elements[i].emap[6]].z +
1373 4 * u * nodes[elements[i].emap[8]].z) +
1374 ((-4 * w * nodes[elements[i].emap[9]].x +
1375 4 * v * (nodes[elements[i].emap[9]].x - nodes[elements[i].emap[2]].x) +
1376 nodes[elements[i].emap[2]].x - nodes[elements[i].emap[3]].x +
1377 4 * w * nodes[elements[i].emap[3]].x -
1378 4 * t * nodes[elements[i].emap[5]].x +
1379 4 * t * nodes[elements[i].emap[6]].x -
1380 4 * u * nodes[elements[i].emap[7]].x +
1381 4 * u * nodes[elements[i].emap[8]].x) *
1382 ((-1 + 4 * t) * nodes[elements[i].emap[0]].y +
1383 4 * (u * nodes[elements[i].emap[4]].y +
1384 v * nodes[elements[i].emap[5]].y +
1385 w * nodes[elements[i].emap[6]].y)) +
1386 ((-1 + 4 * t) * nodes[elements[i].emap[0]].x -
1387 4 * v * nodes[elements[i].emap[9]].x + nodes[elements[i].emap[3]].x -
1388 4 * w * nodes[elements[i].emap[3]].x +
1389 4 * u * nodes[elements[i].emap[4]].x +
1390 4 * v * nodes[elements[i].emap[5]].x -
1391 4 * t * nodes[elements[i].emap[6]].x +
1392 4 * w * nodes[elements[i].emap[6]].x -
1393 4 * u * nodes[elements[i].emap[8]].x) *
1394 (4 * w * nodes[elements[i].emap[9]].y -
1395 nodes[elements[i].emap[2]].y +
1396 4 * v * nodes[elements[i].emap[2]].y +
1397 4 * t * nodes[elements[i].emap[5]].y +
1398 4 * u * nodes[elements[i].emap[7]].y) -
1399 ((-1 + 4 * t) * nodes[elements[i].emap[0]].x -
1400 4 * w * nodes[elements[i].emap[9]].x + nodes[elements[i].emap[2]].x -
1401 4 * v * nodes[elements[i].emap[2]].x +
1402 4 * u * nodes[elements[i].emap[4]].x -
1403 4 * t * nodes[elements[i].emap[5]].x +
1404 4 * v * nodes[elements[i].emap[5]].x +
1405 4 * w * nodes[elements[i].emap[6]].x -
1406 4 * u * nodes[elements[i].emap[7]].x) *
1407 (4 * v * nodes[elements[i].emap[9]].y -
1408 nodes[elements[i].emap[3]].y +
1409 4 * w * nodes[elements[i].emap[3]].y +
1410 4 * t * nodes[elements[i].emap[6]].y +
1411 4 * u * nodes[elements[i].emap[8]].y)) *
1412 ((-1 + 4 * u) * nodes[elements[i].emap[1]].z +
1413 4 * (t * nodes[elements[i].emap[4]].z +
1414 v * nodes[elements[i].emap[7]].z +
1415 w * nodes[elements[i].emap[8]].z));
1416
1417 jac[0][0] =
1418 -((((-1 + 4 * u) * nodes[elements[i].emap[1]].x +
1419 4 * (t * nodes[elements[i].emap[4]].x +
1420 v * nodes[elements[i].emap[7]].x +
1421 w * nodes[elements[i].emap[8]].x)) *
1422 (4 * v * nodes[elements[i].emap[9]].y -
1423 nodes[elements[i].emap[3]].y +
1424 4 * w * nodes[elements[i].emap[3]].y +
1425 4 * t * nodes[elements[i].emap[6]].y +
1426 4 * u * nodes[elements[i].emap[8]].y) -
1427 (4 * v * nodes[elements[i].emap[9]].x - nodes[elements[i].emap[3]].x +
1428 4 * w * nodes[elements[i].emap[3]].x +
1429 4 * t * nodes[elements[i].emap[6]].x +
1430 4 * u * nodes[elements[i].emap[8]].x) *
1431 ((-1 + 4 * u) * nodes[elements[i].emap[1]].y +
1432 4 * (t * nodes[elements[i].emap[4]].y +
1433 v * nodes[elements[i].emap[7]].y +
1434 w * nodes[elements[i].emap[8]].y))) *
1435 (4 * w * nodes[elements[i].emap[9]].z - nodes[elements[i].emap[2]].z +
1436 4 * v * nodes[elements[i].emap[2]].z +
1437 4 * t * nodes[elements[i].emap[5]].z +
1438 4 * u * nodes[elements[i].emap[7]].z)) +
1439 (((-1 + 4 * u) * nodes[elements[i].emap[1]].x +
1440 4 * (t * nodes[elements[i].emap[4]].x +
1441 v * nodes[elements[i].emap[7]].x +
1442 w * nodes[elements[i].emap[8]].x)) *
1443 (4 * w * nodes[elements[i].emap[9]].y -
1444 nodes[elements[i].emap[2]].y +
1445 4 * v * nodes[elements[i].emap[2]].y +
1446 4 * t * nodes[elements[i].emap[5]].y +
1447 4 * u * nodes[elements[i].emap[7]].y) -
1448 (4 * w * nodes[elements[i].emap[9]].x - nodes[elements[i].emap[2]].x +
1449 4 * v * nodes[elements[i].emap[2]].x +
1450 4 * t * nodes[elements[i].emap[5]].x +
1451 4 * u * nodes[elements[i].emap[7]].x) *
1452 ((-1 + 4 * u) * nodes[elements[i].emap[1]].y +
1453 4 * (t * nodes[elements[i].emap[4]].y +
1454 v * nodes[elements[i].emap[7]].y +
1455 w * nodes[elements[i].emap[8]].y))) *
1456 (4 * v * nodes[elements[i].emap[9]].z - nodes[elements[i].emap[3]].z +
1457 4 * w * nodes[elements[i].emap[3]].z +
1458 4 * t * nodes[elements[i].emap[6]].z +
1459 4 * u * nodes[elements[i].emap[8]].z) +
1460 (-((4 * v * nodes[elements[i].emap[9]].x - nodes[elements[i].emap[3]].x +
1461 4 * w * nodes[elements[i].emap[3]].x +
1462 4 * t * nodes[elements[i].emap[6]].x +
1463 4 * u * nodes[elements[i].emap[8]].x) *
1464 (4 * w * nodes[elements[i].emap[9]].y - nodes[elements[i].emap[2]].y +
1465 4 * v * nodes[elements[i].emap[2]].y +
1466 4 * t * nodes[elements[i].emap[5]].y +
1467 4 * u * nodes[elements[i].emap[7]].y)) +
1468 (4 * w * nodes[elements[i].emap[9]].x - nodes[elements[i].emap[2]].x +
1469 4 * v * nodes[elements[i].emap[2]].x +
1470 4 * t * nodes[elements[i].emap[5]].x +
1471 4 * u * nodes[elements[i].emap[7]].x) *
1472 (4 * v * nodes[elements[i].emap[9]].y -
1473 nodes[elements[i].emap[3]].y +
1474 4 * w * nodes[elements[i].emap[3]].y +
1475 4 * t * nodes[elements[i].emap[6]].y +
1476 4 * u * nodes[elements[i].emap[8]].y)) *
1477 ((-1 + 4 * u) * nodes[elements[i].emap[1]].z +
1478 4 * (t * nodes[elements[i].emap[4]].z +
1479 v * nodes[elements[i].emap[7]].z +
1480 w * nodes[elements[i].emap[8]].z));
1481
1482 jac[0][1] =
1483 (nodes[elements[i].emap[1]].y - 4 * u * nodes[elements[i].emap[1]].y -
1484 nodes[elements[i].emap[3]].y + 4 * w * nodes[elements[i].emap[3]].y -
1485 4 * t * nodes[elements[i].emap[4]].y +
1486 4 * t * nodes[elements[i].emap[6]].y +
1487 4 * v * (nodes[elements[i].emap[9]].y - nodes[elements[i].emap[7]].y) +
1488 4 * u * nodes[elements[i].emap[8]].y -
1489 4 * w * nodes[elements[i].emap[8]].y) *
1490 (4 * w * nodes[elements[i].emap[9]].z - nodes[elements[i].emap[2]].z +
1491 4 * v * nodes[elements[i].emap[2]].z +
1492 4 * t * nodes[elements[i].emap[5]].z +
1493 4 * u * nodes[elements[i].emap[7]].z) +
1494 (-4 * w * nodes[elements[i].emap[9]].y - nodes[elements[i].emap[1]].y +
1495 4 * u * nodes[elements[i].emap[1]].y + nodes[elements[i].emap[2]].y -
1496 4 * v * nodes[elements[i].emap[2]].y +
1497 4 * t * nodes[elements[i].emap[4]].y -
1498 4 * t * nodes[elements[i].emap[5]].y -
1499 4 * u * nodes[elements[i].emap[7]].y +
1500 4 * v * nodes[elements[i].emap[7]].y +
1501 4 * w * nodes[elements[i].emap[8]].y) *
1502 (4 * v * nodes[elements[i].emap[9]].z - nodes[elements[i].emap[3]].z +
1503 4 * w * nodes[elements[i].emap[3]].z +
1504 4 * t * nodes[elements[i].emap[6]].z +
1505 4 * u * nodes[elements[i].emap[8]].z) +
1506 (-4 * v * nodes[elements[i].emap[9]].y +
1507 4 * w * nodes[elements[i].emap[9]].y - nodes[elements[i].emap[2]].y +
1508 4 * v * nodes[elements[i].emap[2]].y + nodes[elements[i].emap[3]].y -
1509 4 * w * nodes[elements[i].emap[3]].y +
1510 4 * t * nodes[elements[i].emap[5]].y -
1511 4 * t * nodes[elements[i].emap[6]].y +
1512 4 * u * nodes[elements[i].emap[7]].y -
1513 4 * u * nodes[elements[i].emap[8]].y) *
1514 ((-1 + 4 * u) * nodes[elements[i].emap[1]].z +
1515 4 * (t * nodes[elements[i].emap[4]].z +
1516 v * nodes[elements[i].emap[7]].z +
1517 w * nodes[elements[i].emap[8]].z));
1518
1519 jac[0][2] =
1520 (-4 * v * nodes[elements[i].emap[9]].x - nodes[elements[i].emap[1]].x +
1521 4 * u * nodes[elements[i].emap[1]].x + nodes[elements[i].emap[3]].x -
1522 4 * w * nodes[elements[i].emap[3]].x +
1523 4 * t * nodes[elements[i].emap[4]].x -
1524 4 * t * nodes[elements[i].emap[6]].x +
1525 4 * v * nodes[elements[i].emap[7]].x -
1526 4 * u * nodes[elements[i].emap[8]].x +
1527 4 * w * nodes[elements[i].emap[8]].x) *
1528 (4 * w * nodes[elements[i].emap[9]].z - nodes[elements[i].emap[2]].z +
1529 4 * v * nodes[elements[i].emap[2]].z +
1530 4 * t * nodes[elements[i].emap[5]].z +
1531 4 * u * nodes[elements[i].emap[7]].z) +
1532 (nodes[elements[i].emap[1]].x - 4 * u * nodes[elements[i].emap[1]].x -
1533 nodes[elements[i].emap[2]].x + 4 * v * nodes[elements[i].emap[2]].x -
1534 4 * t * nodes[elements[i].emap[4]].x +
1535 4 * t * nodes[elements[i].emap[5]].x +
1536 4 * u * nodes[elements[i].emap[7]].x -
1537 4 * v * nodes[elements[i].emap[7]].x +
1538 4 * w * (nodes[elements[i].emap[9]].x - nodes[elements[i].emap[8]].x)) *
1539 (4 * v * nodes[elements[i].emap[9]].z - nodes[elements[i].emap[3]].z +
1540 4 * w * nodes[elements[i].emap[3]].z +
1541 4 * t * nodes[elements[i].emap[6]].z +
1542 4 * u * nodes[elements[i].emap[8]].z) +
1543 (-4 * w * nodes[elements[i].emap[9]].x +
1544 4 * v * (nodes[elements[i].emap[9]].x - nodes[elements[i].emap[2]].x) +
1545 nodes[elements[i].emap[2]].x - nodes[elements[i].emap[3]].x +
1546 4 * w * nodes[elements[i].emap[3]].x -
1547 4 * t * nodes[elements[i].emap[5]].x +
1548 4 * t * nodes[elements[i].emap[6]].x -
1549 4 * u * nodes[elements[i].emap[7]].x +
1550 4 * u * nodes[elements[i].emap[8]].x) *
1551 ((-1 + 4 * u) * nodes[elements[i].emap[1]].z +
1552 4 * (t * nodes[elements[i].emap[4]].z +
1553 v * nodes[elements[i].emap[7]].z +
1554 w * nodes[elements[i].emap[8]].z));
1555
1556 jac[0][3] =
1557 (nodes[elements[i].emap[1]].x - 4 * u * nodes[elements[i].emap[1]].x -
1558 nodes[elements[i].emap[3]].x + 4 * w * nodes[elements[i].emap[3]].x -
1559 4 * t * nodes[elements[i].emap[4]].x +
1560 4 * t * nodes[elements[i].emap[6]].x +
1561 4 * v * (nodes[elements[i].emap[9]].x - nodes[elements[i].emap[7]].x) +
1562 4 * u * nodes[elements[i].emap[8]].x -
1563 4 * w * nodes[elements[i].emap[8]].x) *
1564 (4 * w * nodes[elements[i].emap[9]].y - nodes[elements[i].emap[2]].y +
1565 4 * v * nodes[elements[i].emap[2]].y +
1566 4 * t * nodes[elements[i].emap[5]].y +
1567 4 * u * nodes[elements[i].emap[7]].y) +
1568 (-4 * w * nodes[elements[i].emap[9]].x - nodes[elements[i].emap[1]].x +
1569 4 * u * nodes[elements[i].emap[1]].x + nodes[elements[i].emap[2]].x -
1570 4 * v * nodes[elements[i].emap[2]].x +
1571 4 * t * nodes[elements[i].emap[4]].x -
1572 4 * t * nodes[elements[i].emap[5]].x -
1573 4 * u * nodes[elements[i].emap[7]].x +
1574 4 * v * nodes[elements[i].emap[7]].x +
1575 4 * w * nodes[elements[i].emap[8]].x) *
1576 (4 * v * nodes[elements[i].emap[9]].y - nodes[elements[i].emap[3]].y +
1577 4 * w * nodes[elements[i].emap[3]].y +
1578 4 * t * nodes[elements[i].emap[6]].y +
1579 4 * u * nodes[elements[i].emap[8]].y) +
1580 (-4 * v * nodes[elements[i].emap[9]].x +
1581 4 * w * nodes[elements[i].emap[9]].x - nodes[elements[i].emap[2]].x +
1582 4 * v * nodes[elements[i].emap[2]].x + nodes[elements[i].emap[3]].x -
1583 4 * w * nodes[elements[i].emap[3]].x +
1584 4 * t * nodes[elements[i].emap[5]].x -
1585 4 * t * nodes[elements[i].emap[6]].x +
1586 4 * u * nodes[elements[i].emap[7]].x -
1587 4 * u * nodes[elements[i].emap[8]].x) *
1588 ((-1 + 4 * u) * nodes[elements[i].emap[1]].y +
1589 4 * (t * nodes[elements[i].emap[4]].y +
1590 v * nodes[elements[i].emap[7]].y +
1591 w * nodes[elements[i].emap[8]].y));
1592
1593 jac[1][0] =
1594 -((-((4 * v * nodes[elements[i].emap[9]].x -
1595 nodes[elements[i].emap[3]].x +
1596 4 * w * nodes[elements[i].emap[3]].x +
1597 4 * t * nodes[elements[i].emap[6]].x +
1598 4 * u * nodes[elements[i].emap[8]].x) *
1599 (4 * w * nodes[elements[i].emap[9]].y -
1600 nodes[elements[i].emap[2]].y +
1601 4 * v * nodes[elements[i].emap[2]].y +
1602 4 * t * nodes[elements[i].emap[5]].y +
1603 4 * u * nodes[elements[i].emap[7]].y)) +
1604 (4 * w * nodes[elements[i].emap[9]].x - nodes[elements[i].emap[2]].x +
1605 4 * v * nodes[elements[i].emap[2]].x +
1606 4 * t * nodes[elements[i].emap[5]].x +
1607 4 * u * nodes[elements[i].emap[7]].x) *
1608 (4 * v * nodes[elements[i].emap[9]].y -
1609 nodes[elements[i].emap[3]].y +
1610 4 * w * nodes[elements[i].emap[3]].y +
1611 4 * t * nodes[elements[i].emap[6]].y +
1612 4 * u * nodes[elements[i].emap[8]].y)) *
1613 ((-1 + 4 * t) * nodes[elements[i].emap[0]].z +
1614 4 * (u * nodes[elements[i].emap[4]].z +
1615 v * nodes[elements[i].emap[5]].z +
1616 w * nodes[elements[i].emap[6]].z))) +
1617 (-((4 * v * nodes[elements[i].emap[9]].x - nodes[elements[i].emap[3]].x +
1618 4 * w * nodes[elements[i].emap[3]].x +
1619 4 * t * nodes[elements[i].emap[6]].x +
1620 4 * u * nodes[elements[i].emap[8]].x) *
1621 ((-1 + 4 * t) * nodes[elements[i].emap[0]].y +
1622 4 * (u * nodes[elements[i].emap[4]].y +
1623 v * nodes[elements[i].emap[5]].y +
1624 w * nodes[elements[i].emap[6]].y))) +
1625 ((-1 + 4 * t) * nodes[elements[i].emap[0]].x +
1626 4 * (u * nodes[elements[i].emap[4]].x +
1627 v * nodes[elements[i].emap[5]].x +
1628 w * nodes[elements[i].emap[6]].x)) *
1629 (4 * v * nodes[elements[i].emap[9]].y -
1630 nodes[elements[i].emap[3]].y +
1631 4 * w * nodes[elements[i].emap[3]].y +
1632 4 * t * nodes[elements[i].emap[6]].y +
1633 4 * u * nodes[elements[i].emap[8]].y)) *
1634 (4 * w * nodes[elements[i].emap[9]].z - nodes[elements[i].emap[2]].z +
1635 4 * v * nodes[elements[i].emap[2]].z +
1636 4 * t * nodes[elements[i].emap[5]].z +
1637 4 * u * nodes[elements[i].emap[7]].z) -
1638 (-((4 * w * nodes[elements[i].emap[9]].x - nodes[elements[i].emap[2]].x +
1639 4 * v * nodes[elements[i].emap[2]].x +
1640 4 * t * nodes[elements[i].emap[5]].x +
1641 4 * u * nodes[elements[i].emap[7]].x) *
1642 ((-1 + 4 * t) * nodes[elements[i].emap[0]].y +
1643 4 * (u * nodes[elements[i].emap[4]].y +
1644 v * nodes[elements[i].emap[5]].y +
1645 w * nodes[elements[i].emap[6]].y))) +
1646 ((-1 + 4 * t) * nodes[elements[i].emap[0]].x +
1647 4 * (u * nodes[elements[i].emap[4]].x +
1648 v * nodes[elements[i].emap[5]].x +
1649 w * nodes[elements[i].emap[6]].x)) *
1650 (4 * w * nodes[elements[i].emap[9]].y -
1651 nodes[elements[i].emap[2]].y +
1652 4 * v * nodes[elements[i].emap[2]].y +
1653 4 * t * nodes[elements[i].emap[5]].y +
1654 4 * u * nodes[elements[i].emap[7]].y)) *
1655 (4 * v * nodes[elements[i].emap[9]].z - nodes[elements[i].emap[3]].z +
1656 4 * w * nodes[elements[i].emap[3]].z +
1657 4 * t * nodes[elements[i].emap[6]].z +
1658 4 * u * nodes[elements[i].emap[8]].z);
1659
1660 jac[1][1] =
1661 (-4 * w * nodes[elements[i].emap[9]].y +
1662 4 * v * (nodes[elements[i].emap[9]].y - nodes[elements[i].emap[2]].y) +
1663 nodes[elements[i].emap[2]].y - nodes[elements[i].emap[3]].y +
1664 4 * w * nodes[elements[i].emap[3]].y -
1665 4 * t * nodes[elements[i].emap[5]].y +
1666 4 * t * nodes[elements[i].emap[6]].y -
1667 4 * u * nodes[elements[i].emap[7]].y +
1668 4 * u * nodes[elements[i].emap[8]].y) *
1669 ((-1 + 4 * t) * nodes[elements[i].emap[0]].z +
1670 4 * (u * nodes[elements[i].emap[4]].z +
1671 v * nodes[elements[i].emap[5]].z +
1672 w * nodes[elements[i].emap[6]].z)) +
1673 ((-1 + 4 * t) * nodes[elements[i].emap[0]].y -
1674 4 * v * nodes[elements[i].emap[9]].y + nodes[elements[i].emap[3]].y -
1675 4 * w * nodes[elements[i].emap[3]].y +
1676 4 * u * nodes[elements[i].emap[4]].y +
1677 4 * v * nodes[elements[i].emap[5]].y -
1678 4 * t * nodes[elements[i].emap[6]].y +
1679 4 * w * nodes[elements[i].emap[6]].y -
1680 4 * u * nodes[elements[i].emap[8]].y) *
1681 (4 * w * nodes[elements[i].emap[9]].z - nodes[elements[i].emap[2]].z +
1682 4 * v * nodes[elements[i].emap[2]].z +
1683 4 * t * nodes[elements[i].emap[5]].z +
1684 4 * u * nodes[elements[i].emap[7]].z) -
1685 ((-1 + 4 * t) * nodes[elements[i].emap[0]].y -
1686 4 * w * nodes[elements[i].emap[9]].y + nodes[elements[i].emap[2]].y -
1687 4 * v * nodes[elements[i].emap[2]].y +
1688 4 * u * nodes[elements[i].emap[4]].y -
1689 4 * t * nodes[elements[i].emap[5]].y +
1690 4 * v * nodes[elements[i].emap[5]].y +
1691 4 * w * nodes[elements[i].emap[6]].y -
1692 4 * u * nodes[elements[i].emap[7]].y) *
1693 (4 * v * nodes[elements[i].emap[9]].z - nodes[elements[i].emap[3]].z +
1694 4 * w * nodes[elements[i].emap[3]].z +
1695 4 * t * nodes[elements[i].emap[6]].z +
1696 4 * u * nodes[elements[i].emap[8]].z);
1697
1698 jac[1][2] =
1699 (-4 * v * nodes[elements[i].emap[9]].x +
1700 4 * w * nodes[elements[i].emap[9]].x - nodes[elements[i].emap[2]].x +
1701 4 * v * nodes[elements[i].emap[2]].x + nodes[elements[i].emap[3]].x -
1702 4 * w * nodes[elements[i].emap[3]].x +
1703 4 * t * nodes[elements[i].emap[5]].x -
1704 4 * t * nodes[elements[i].emap[6]].x +
1705 4 * u * nodes[elements[i].emap[7]].x -
1706 4 * u * nodes[elements[i].emap[8]].x) *
1707 ((-1 + 4 * t) * nodes[elements[i].emap[0]].z +
1708 4 * (u * nodes[elements[i].emap[4]].z +
1709 v * nodes[elements[i].emap[5]].z +
1710 w * nodes[elements[i].emap[6]].z)) -
1711 ((-1 + 4 * t) * nodes[elements[i].emap[0]].x -
1712 4 * v * nodes[elements[i].emap[9]].x + nodes[elements[i].emap[3]].x -
1713 4 * w * nodes[elements[i].emap[3]].x +
1714 4 * u * nodes[elements[i].emap[4]].x +
1715 4 * v * nodes[elements[i].emap[5]].x -
1716 4 * t * nodes[elements[i].emap[6]].x +
1717 4 * w * nodes[elements[i].emap[6]].x -
1718 4 * u * nodes[elements[i].emap[8]].x) *
1719 (4 * w * nodes[elements[i].emap[9]].z - nodes[elements[i].emap[2]].z +
1720 4 * v * nodes[elements[i].emap[2]].z +
1721 4 * t * nodes[elements[i].emap[5]].z +
1722 4 * u * nodes[elements[i].emap[7]].z) +
1723 ((-1 + 4 * t) * nodes[elements[i].emap[0]].x -
1724 4 * w * nodes[elements[i].emap[9]].x + nodes[elements[i].emap[2]].x -
1725 4 * v * nodes[elements[i].emap[2]].x +
1726 4 * u * nodes[elements[i].emap[4]].x -
1727 4 * t * nodes[elements[i].emap[5]].x +
1728 4 * v * nodes[elements[i].emap[5]].x +
1729 4 * w * nodes[elements[i].emap[6]].x -
1730 4 * u * nodes[elements[i].emap[7]].x) *
1731 (4 * v * nodes[elements[i].emap[9]].z - nodes[elements[i].emap[3]].z +
1732 4 * w * nodes[elements[i].emap[3]].z +
1733 4 * t * nodes[elements[i].emap[6]].z +
1734 4 * u * nodes[elements[i].emap[8]].z);
1735
1736 jac[1][3] =
1737 (-4 * w * nodes[elements[i].emap[9]].x +
1738 4 * v * (nodes[elements[i].emap[9]].x - nodes[elements[i].emap[2]].x) +
1739 nodes[elements[i].emap[2]].x - nodes[elements[i].emap[3]].x +
1740 4 * w * nodes[elements[i].emap[3]].x -
1741 4 * t * nodes[elements[i].emap[5]].x +
1742 4 * t * nodes[elements[i].emap[6]].x -
1743 4 * u * nodes[elements[i].emap[7]].x +
1744 4 * u * nodes[elements[i].emap[8]].x) *
1745 ((-1 + 4 * t) * nodes[elements[i].emap[0]].y +
1746 4 * (u * nodes[elements[i].emap[4]].y +
1747 v * nodes[elements[i].emap[5]].y +
1748 w * nodes[elements[i].emap[6]].y)) +
1749 ((-1 + 4 * t) * nodes[elements[i].emap[0]].x -
1750 4 * v * nodes[elements[i].emap[9]].x + nodes[elements[i].emap[3]].x -
1751 4 * w * nodes[elements[i].emap[3]].x +
1752 4 * u * nodes[elements[i].emap[4]].x +
1753 4 * v * nodes[elements[i].emap[5]].x -
1754 4 * t * nodes[elements[i].emap[6]].x +
1755 4 * w * nodes[elements[i].emap[6]].x -
1756 4 * u * nodes[elements[i].emap[8]].x) *
1757 (4 * w * nodes[elements[i].emap[9]].y - nodes[elements[i].emap[2]].y +
1758 4 * v * nodes[elements[i].emap[2]].y +
1759 4 * t * nodes[elements[i].emap[5]].y +
1760 4 * u * nodes[elements[i].emap[7]].y) -
1761 ((-1 + 4 * t) * nodes[elements[i].emap[0]].x -
1762 4 * w * nodes[elements[i].emap[9]].x + nodes[elements[i].emap[2]].x -
1763 4 * v * nodes[elements[i].emap[2]].x +
1764 4 * u * nodes[elements[i].emap[4]].x -
1765 4 * t * nodes[elements[i].emap[5]].x +
1766 4 * v * nodes[elements[i].emap[5]].x +
1767 4 * w * nodes[elements[i].emap[6]].x -
1768 4 * u * nodes[elements[i].emap[7]].x) *
1769 (4 * v * nodes[elements[i].emap[9]].y - nodes[elements[i].emap[3]].y +
1770 4 * w * nodes[elements[i].emap[3]].y +
1771 4 * t * nodes[elements[i].emap[6]].y +
1772 4 * u * nodes[elements[i].emap[8]].y);
1773
1774 jac[2][0] =
1775 (((-1 + 4 * u) * nodes[elements[i].emap[1]].x +
1776 4 * (t * nodes[elements[i].emap[4]].x +
1777 v * nodes[elements[i].emap[7]].x +
1778 w * nodes[elements[i].emap[8]].x)) *
1779 (4 * v * nodes[elements[i].emap[9]].y -
1780 nodes[elements[i].emap[3]].y +
1781 4 * w * nodes[elements[i].emap[3]].y +
1782 4 * t * nodes[elements[i].emap[6]].y +
1783 4 * u * nodes[elements[i].emap[8]].y) -
1784 (4 * v * nodes[elements[i].emap[9]].x - nodes[elements[i].emap[3]].x +
1785 4 * w * nodes[elements[i].emap[3]].x +
1786 4 * t * nodes[elements[i].emap[6]].x +
1787 4 * u * nodes[elements[i].emap[8]].x) *
1788 ((-1 + 4 * u) * nodes[elements[i].emap[1]].y +
1789 4 * (t * nodes[elements[i].emap[4]].y +
1790 v * nodes[elements[i].emap[7]].y +
1791 w * nodes[elements[i].emap[8]].y))) *
1792 ((-1 + 4 * t) * nodes[elements[i].emap[0]].z +
1793 4 * (u * nodes[elements[i].emap[4]].z +
1794 v * nodes[elements[i].emap[5]].z +
1795 w * nodes[elements[i].emap[6]].z)) +
1796 (-(((-1 + 4 * u) * nodes[elements[i].emap[1]].x +
1797 4 * (t * nodes[elements[i].emap[4]].x +
1798 v * nodes[elements[i].emap[7]].x +
1799 w * nodes[elements[i].emap[8]].x)) *
1800 ((-1 + 4 * t) * nodes[elements[i].emap[0]].y +
1801 4 * (u * nodes[elements[i].emap[4]].y +
1802 v * nodes[elements[i].emap[5]].y +
1803 w * nodes[elements[i].emap[6]].y))) +
1804 ((-1 + 4 * t) * nodes[elements[i].emap[0]].x +
1805 4 * (u * nodes[elements[i].emap[4]].x +
1806 v * nodes[elements[i].emap[5]].x +
1807 w * nodes[elements[i].emap[6]].x)) *
1808 ((-1 + 4 * u) * nodes[elements[i].emap[1]].y +
1809 4 * (t * nodes[elements[i].emap[4]].y +
1810 v * nodes[elements[i].emap[7]].y +
1811 w * nodes[elements[i].emap[8]].y))) *
1812 (4 * v * nodes[elements[i].emap[9]].z - nodes[elements[i].emap[3]].z +
1813 4 * w * nodes[elements[i].emap[3]].z +
1814 4 * t * nodes[elements[i].emap[6]].z +
1815 4 * u * nodes[elements[i].emap[8]].z) -
1816 (-((4 * v * nodes[elements[i].emap[9]].x - nodes[elements[i].emap[3]].x +
1817 4 * w * nodes[elements[i].emap[3]].x +
1818 4 * t * nodes[elements[i].emap[6]].x +
1819 4 * u * nodes[elements[i].emap[8]].x) *
1820 ((-1 + 4 * t) * nodes[elements[i].emap[0]].y +
1821 4 * (u * nodes[elements[i].emap[4]].y +
1822 v * nodes[elements[i].emap[5]].y +
1823 w * nodes[elements[i].emap[6]].y))) +
1824 ((-1 + 4 * t) * nodes[elements[i].emap[0]].x +
1825 4 * (u * nodes[elements[i].emap[4]].x +
1826 v * nodes[elements[i].emap[5]].x +
1827 w * nodes[elements[i].emap[6]].x)) *
1828 (4 * v * nodes[elements[i].emap[9]].y -
1829 nodes[elements[i].emap[3]].y +
1830 4 * w * nodes[elements[i].emap[3]].y +
1831 4 * t * nodes[elements[i].emap[6]].y +
1832 4 * u * nodes[elements[i].emap[8]].y)) *
1833 ((-1 + 4 * u) * nodes[elements[i].emap[1]].z +
1834 4 * (t * nodes[elements[i].emap[4]].z +
1835 v * nodes[elements[i].emap[7]].z +
1836 w * nodes[elements[i].emap[8]].z));
1837
1838 jac[2][1] =
1839 (-4 * v * nodes[elements[i].emap[9]].y - nodes[elements[i].emap[1]].y +
1840 4 * u * nodes[elements[i].emap[1]].y + nodes[elements[i].emap[3]].y -
1841 4 * w * nodes[elements[i].emap[3]].y +
1842 4 * t * nodes[elements[i].emap[4]].y -
1843 4 * t * nodes[elements[i].emap[6]].y +
1844 4 * v * nodes[elements[i].emap[7]].y -
1845 4 * u * nodes[elements[i].emap[8]].y +
1846 4 * w * nodes[elements[i].emap[8]].y) *
1847 ((-1 + 4 * t) * nodes[elements[i].emap[0]].z +
1848 4 * (u * nodes[elements[i].emap[4]].z +
1849 v * nodes[elements[i].emap[5]].z +
1850 w * nodes[elements[i].emap[6]].z)) +
1851 ((-1 + 4 * t) * nodes[elements[i].emap[0]].y +
1852 nodes[elements[i].emap[1]].y - 4 * u * nodes[elements[i].emap[1]].y +
1853 4 * (-(t * nodes[elements[i].emap[4]].y) +
1854 u * nodes[elements[i].emap[4]].y +
1855 v * nodes[elements[i].emap[5]].y +
1856 w * nodes[elements[i].emap[6]].y -
1857 v * nodes[elements[i].emap[7]].y -
1858 w * nodes[elements[i].emap[8]].y)) *
1859 (4 * v * nodes[elements[i].emap[9]].z - nodes[elements[i].emap[3]].z +
1860 4 * w * nodes[elements[i].emap[3]].z +
1861 4 * t * nodes[elements[i].emap[6]].z +
1862 4 * u * nodes[elements[i].emap[8]].z) -
1863 ((-1 + 4 * t) * nodes[elements[i].emap[0]].y -
1864 4 * v * nodes[elements[i].emap[9]].y + nodes[elements[i].emap[3]].y -
1865 4 * w * nodes[elements[i].emap[3]].y +
1866 4 * u * nodes[elements[i].emap[4]].y +
1867 4 * v * nodes[elements[i].emap[5]].y -
1868 4 * t * nodes[elements[i].emap[6]].y +
1869 4 * w * nodes[elements[i].emap[6]].y -
1870 4 * u * nodes[elements[i].emap[8]].y) *
1871 ((-1 + 4 * u) * nodes[elements[i].emap[1]].z +
1872 4 * (t * nodes[elements[i].emap[4]].z +
1873 v * nodes[elements[i].emap[7]].z +
1874 w * nodes[elements[i].emap[8]].z));
1875
1876 jac[2][2] =
1877 (nodes[elements[i].emap[1]].x - 4 * u * nodes[elements[i].emap[1]].x -
1878 nodes[elements[i].emap[3]].x + 4 * w * nodes[elements[i].emap[3]].x -
1879 4 * t * nodes[elements[i].emap[4]].x +
1880 4 * t * nodes[elements[i].emap[6]].x +
1881 4 * v * (nodes[elements[i].emap[9]].x - nodes[elements[i].emap[7]].x) +
1882 4 * u * nodes[elements[i].emap[8]].x -
1883 4 * w * nodes[elements[i].emap[8]].x) *
1884 ((-1 + 4 * t) * nodes[elements[i].emap[0]].z +
1885 4 * (u * nodes[elements[i].emap[4]].z +
1886 v * nodes[elements[i].emap[5]].z +
1887 w * nodes[elements[i].emap[6]].z)) -
1888 ((-1 + 4 * t) * nodes[elements[i].emap[0]].x +
1889 nodes[elements[i].emap[1]].x - 4 * u * nodes[elements[i].emap[1]].x +
1890 4 * (-(t * nodes[elements[i].emap[4]].x) +
1891 u * nodes[elements[i].emap[4]].x +
1892 v * nodes[elements[i].emap[5]].x +
1893 w * nodes[elements[i].emap[6]].x -
1894 v * nodes[elements[i].emap[7]].x -
1895 w * nodes[elements[i].emap[8]].x)) *
1896 (4 * v * nodes[elements[i].emap[9]].z - nodes[elements[i].emap[3]].z +
1897 4 * w * nodes[elements[i].emap[3]].z +
1898 4 * t * nodes[elements[i].emap[6]].z +
1899 4 * u * nodes[elements[i].emap[8]].z) +
1900 ((-1 + 4 * t) * nodes[elements[i].emap[0]].x -
1901 4 * v * nodes[elements[i].emap[9]].x + nodes[elements[i].emap[3]].x -
1902 4 * w * nodes[elements[i].emap[3]].x +
1903 4 * u * nodes[elements[i].emap[4]].x +
1904 4 * v * nodes[elements[i].emap[5]].x -
1905 4 * t * nodes[elements[i].emap[6]].x +
1906 4 * w * nodes[elements[i].emap[6]].x -
1907 4 * u * nodes[elements[i].emap[8]].x) *
1908 ((-1 + 4 * u) * nodes[elements[i].emap[1]].z +
1909 4 * (t * nodes[elements[i].emap[4]].z +
1910 v * nodes[elements[i].emap[7]].z +
1911 w * nodes[elements[i].emap[8]].z));
1912
1913 jac[2][3] =
1914 (-4 * v * nodes[elements[i].emap[9]].x - nodes[elements[i].emap[1]].x +
1915 4 * u * nodes[elements[i].emap[1]].x + nodes[elements[i].emap[3]].x -
1916 4 * w * nodes[elements[i].emap[3]].x +
1917 4 * t * nodes[elements[i].emap[4]].x -
1918 4 * t * nodes[elements[i].emap[6]].x +
1919 4 * v * nodes[elements[i].emap[7]].x -
1920 4 * u * nodes[elements[i].emap[8]].x +
1921 4 * w * nodes[elements[i].emap[8]].x) *
1922 ((-1 + 4 * t) * nodes[elements[i].emap[0]].y +
1923 4 * (u * nodes[elements[i].emap[4]].y +
1924 v * nodes[elements[i].emap[5]].y +
1925 w * nodes[elements[i].emap[6]].y)) +
1926 ((-1 + 4 * t) * nodes[elements[i].emap[0]].x +
1927 nodes[elements[i].emap[1]].x - 4 * u * nodes[elements[i].emap[1]].x +
1928 4 * (-(t * nodes[elements[i].emap[4]].x) +
1929 u * nodes[elements[i].emap[4]].x +
1930 v * nodes[elements[i].emap[5]].x +
1931 w * nodes[elements[i].emap[6]].x -
1932 v * nodes[elements[i].emap[7]].x -
1933 w * nodes[elements[i].emap[8]].x)) *
1934 (4 * v * nodes[elements[i].emap[9]].y - nodes[elements[i].emap[3]].y +
1935 4 * w * nodes[elements[i].emap[3]].y +
1936 4 * t * nodes[elements[i].emap[6]].y +
1937 4 * u * nodes[elements[i].emap[8]].y) -
1938 ((-1 + 4 * t) * nodes[elements[i].emap[0]].x -
1939 4 * v * nodes[elements[i].emap[9]].x + nodes[elements[i].emap[3]].x -
1940 4 * w * nodes[elements[i].emap[3]].x +
1941 4 * u * nodes[elements[i].emap[4]].x +
1942 4 * v * nodes[elements[i].emap[5]].x -
1943 4 * t * nodes[elements[i].emap[6]].x +
1944 4 * w * nodes[elements[i].emap[6]].x -
1945 4 * u * nodes[elements[i].emap[8]].x) *
1946 ((-1 + 4 * u) * nodes[elements[i].emap[1]].y +
1947 4 * (t * nodes[elements[i].emap[4]].y +
1948 v * nodes[elements[i].emap[7]].y +
1949 w * nodes[elements[i].emap[8]].y));
1950
1951 jac[3][0] =
1952 -((((-1 + 4 * u) * nodes[elements[i].emap[1]].x +
1953 4 * (t * nodes[elements[i].emap[4]].x +
1954 v * nodes[elements[i].emap[7]].x +
1955 w * nodes[elements[i].emap[8]].x)) *
1956 (4 * w * nodes[elements[i].emap[9]].y -
1957 nodes[elements[i].emap[2]].y +
1958 4 * v * nodes[elements[i].emap[2]].y +
1959 4 * t * nodes[elements[i].emap[5]].y +
1960 4 * u * nodes[elements[i].emap[7]].y) -
1961 (4 * w * nodes[elements[i].emap[9]].x - nodes[elements[i].emap[2]].x +
1962 4 * v * nodes[elements[i].emap[2]].x +
1963 4 * t * nodes[elements[i].emap[5]].x +
1964 4 * u * nodes[elements[i].emap[7]].x) *
1965 ((-1 + 4 * u) * nodes[elements[i].emap[1]].y +
1966 4 * (t * nodes[elements[i].emap[4]].y +
1967 v * nodes[elements[i].emap[7]].y +
1968 w * nodes[elements[i].emap[8]].y))) *
1969 ((-1 + 4 * t) * nodes[elements[i].emap[0]].z +
1970 4 * (u * nodes[elements[i].emap[4]].z +
1971 v * nodes[elements[i].emap[5]].z +
1972 w * nodes[elements[i].emap[6]].z))) -
1973 (-(((-1 + 4 * u) * nodes[elements[i].emap[1]].x +
1974 4 * (t * nodes[elements[i].emap[4]].x +
1975 v * nodes[elements[i].emap[7]].x +
1976 w * nodes[elements[i].emap[8]].x)) *
1977 ((-1 + 4 * t) * nodes[elements[i].emap[0]].y +
1978 4 * (u * nodes[elements[i].emap[4]].y +
1979 v * nodes[elements[i].emap[5]].y +
1980 w * nodes[elements[i].emap[6]].y))) +
1981 ((-1 + 4 * t) * nodes[elements[i].emap[0]].x +
1982 4 * (u * nodes[elements[i].emap[4]].x +
1983 v * nodes[elements[i].emap[5]].x +
1984 w * nodes[elements[i].emap[6]].x)) *
1985 ((-1 + 4 * u) * nodes[elements[i].emap[1]].y +
1986 4 * (t * nodes[elements[i].emap[4]].y +
1987 v * nodes[elements[i].emap[7]].y +
1988 w * nodes[elements[i].emap[8]].y))) *
1989 (4 * w * nodes[elements[i].emap[9]].z - nodes[elements[i].emap[2]].z +
1990 4 * v * nodes[elements[i].emap[2]].z +
1991 4 * t * nodes[elements[i].emap[5]].z +
1992 4 * u * nodes[elements[i].emap[7]].z) +
1993 (-((4 * w * nodes[elements[i].emap[9]].x - nodes[elements[i].emap[2]].x +
1994 4 * v * nodes[elements[i].emap[2]].x +
1995 4 * t * nodes[elements[i].emap[5]].x +
1996 4 * u * nodes[elements[i].emap[7]].x) *
1997 ((-1 + 4 * t) * nodes[elements[i].emap[0]].y +
1998 4 * (u * nodes[elements[i].emap[4]].y +
1999 v * nodes[elements[i].emap[5]].y +
2000 w * nodes[elements[i].emap[6]].y))) +
2001 ((-1 + 4 * t) * nodes[elements[i].emap[0]].x +
2002 4 * (u * nodes[elements[i].emap[4]].x +
2003 v * nodes[elements[i].emap[5]].x +
2004 w * nodes[elements[i].emap[6]].x)) *
2005 (4 * w * nodes[elements[i].emap[9]].y -
2006 nodes[elements[i].emap[2]].y +
2007 4 * v * nodes[elements[i].emap[2]].y +
2008 4 * t * nodes[elements[i].emap[5]].y +
2009 4 * u * nodes[elements[i].emap[7]].y)) *
2010 ((-1 + 4 * u) * nodes[elements[i].emap[1]].z +
2011 4 * (t * nodes[elements[i].emap[4]].z +
2012 v * nodes[elements[i].emap[7]].z +
2013 w * nodes[elements[i].emap[8]].z));
2014
2015 jac[3][1] =
2016 (nodes[elements[i].emap[1]].y - 4 * u * nodes[elements[i].emap[1]].y -
2017 nodes[elements[i].emap[2]].y + 4 * v * nodes[elements[i].emap[2]].y -
2018 4 * t * nodes[elements[i].emap[4]].y +
2019 4 * t * nodes[elements[i].emap[5]].y +
2020 4 * u * nodes[elements[i].emap[7]].y -
2021 4 * v * nodes[elements[i].emap[7]].y +
2022 4 * w * (nodes[elements[i].emap[9]].y - nodes[elements[i].emap[8]].y)) *
2023 ((-1 + 4 * t) * nodes[elements[i].emap[0]].z +
2024 4 * (u * nodes[elements[i].emap[4]].z +
2025 v * nodes[elements[i].emap[5]].z +
2026 w * nodes[elements[i].emap[6]].z)) -
2027 ((-1 + 4 * t) * nodes[elements[i].emap[0]].y +
2028 nodes[elements[i].emap[1]].y - 4 * u * nodes[elements[i].emap[1]].y +
2029 4 * (-(t * nodes[elements[i].emap[4]].y) +
2030 u * nodes[elements[i].emap[4]].y +
2031 v * nodes[elements[i].emap[5]].y +
2032 w * nodes[elements[i].emap[6]].y -
2033 v * nodes[elements[i].emap[7]].y -
2034 w * nodes[elements[i].emap[8]].y)) *
2035 (4 * w * nodes[elements[i].emap[9]].z - nodes[elements[i].emap[2]].z +
2036 4 * v * nodes[elements[i].emap[2]].z +
2037 4 * t * nodes[elements[i].emap[5]].z +
2038 4 * u * nodes[elements[i].emap[7]].z) +
2039 ((-1 + 4 * t) * nodes[elements[i].emap[0]].y -
2040 4 * w * nodes[elements[i].emap[9]].y + nodes[elements[i].emap[2]].y -
2041 4 * v * nodes[elements[i].emap[2]].y +
2042 4 * u * nodes[elements[i].emap[4]].y -
2043 4 * t * nodes[elements[i].emap[5]].y +
2044 4 * v * nodes[elements[i].emap[5]].y +
2045 4 * w * nodes[elements[i].emap[6]].y -
2046 4 * u * nodes[elements[i].emap[7]].y) *
2047 ((-1 + 4 * u) * nodes[elements[i].emap[1]].z +
2048 4 * (t * nodes[elements[i].emap[4]].z +
2049 v * nodes[elements[i].emap[7]].z +
2050 w * nodes[elements[i].emap[8]].z));
2051
2052 jac[3][2] =
2053 (-4 * w * nodes[elements[i].emap[9]].x - nodes[elements[i].emap[1]].x +
2054 4 * u * nodes[elements[i].emap[1]].x + nodes[elements[i].emap[2]].x -
2055 4 * v * nodes[elements[i].emap[2]].x +
2056 4 * t * nodes[elements[i].emap[4]].x -
2057 4 * t * nodes[elements[i].emap[5]].x -
2058 4 * u * nodes[elements[i].emap[7]].x +
2059 4 * v * nodes[elements[i].emap[7]].x +
2060 4 * w * nodes[elements[i].emap[8]].x) *
2061 ((-1 + 4 * t) * nodes[elements[i].emap[0]].z +
2062 4 * (u * nodes[elements[i].emap[4]].z +
2063 v * nodes[elements[i].emap[5]].z +
2064 w * nodes[elements[i].emap[6]].z)) +
2065 ((-1 + 4 * t) * nodes[elements[i].emap[0]].x +
2066 nodes[elements[i].emap[1]].x - 4 * u * nodes[elements[i].emap[1]].x +
2067 4 * (-(t * nodes[elements[i].emap[4]].x) +
2068 u * nodes[elements[i].emap[4]].x +
2069 v * nodes[elements[i].emap[5]].x +
2070 w * nodes[elements[i].emap[6]].x -
2071 v * nodes[elements[i].emap[7]].x -
2072 w * nodes[elements[i].emap[8]].x)) *
2073 (4 * w * nodes[elements[i].emap[9]].z - nodes[elements[i].emap[2]].z +
2074 4 * v * nodes[elements[i].emap[2]].z +
2075 4 * t * nodes[elements[i].emap[5]].z +
2076 4 * u * nodes[elements[i].emap[7]].z) -
2077 ((-1 + 4 * t) * nodes[elements[i].emap[0]].x -
2078 4 * w * nodes[elements[i].emap[9]].x + nodes[elements[i].emap[2]].x -
2079 4 * v * nodes[elements[i].emap[2]].x +
2080 4 * u * nodes[elements[i].emap[4]].x -
2081 4 * t * nodes[elements[i].emap[5]].x +
2082 4 * v * nodes[elements[i].emap[5]].x +
2083 4 * w * nodes[elements[i].emap[6]].x -
2084 4 * u * nodes[elements[i].emap[7]].x) *
2085 ((-1 + 4 * u) * nodes[elements[i].emap[1]].z +
2086 4 * (t * nodes[elements[i].emap[4]].z +
2087 v * nodes[elements[i].emap[7]].z +
2088 w * nodes[elements[i].emap[8]].z));
2089
2090 jac[3][3] =
2091 (nodes[elements[i].emap[1]].x - 4 * u * nodes[elements[i].emap[1]].x -
2092 nodes[elements[i].emap[2]].x + 4 * v * nodes[elements[i].emap[2]].x -
2093 4 * t * nodes[elements[i].emap[4]].x +
2094 4 * t * nodes[elements[i].emap[5]].x +
2095 4 * u * nodes[elements[i].emap[7]].x -
2096 4 * v * nodes[elements[i].emap[7]].x +
2097 4 * w * (nodes[elements[i].emap[9]].x - nodes[elements[i].emap[8]].x)) *
2098 ((-1 + 4 * t) * nodes[elements[i].emap[0]].y +
2099 4 * (u * nodes[elements[i].emap[4]].y +
2100 v * nodes[elements[i].emap[5]].y +
2101 w * nodes[elements[i].emap[6]].y)) -
2102 ((-1 + 4 * t) * nodes[elements[i].emap[0]].x +
2103 nodes[elements[i].emap[1]].x - 4 * u * nodes[elements[i].emap[1]].x +
2104 4 * (-(t * nodes[elements[i].emap[4]].x) +
2105 u * nodes[elements[i].emap[4]].x +
2106 v * nodes[elements[i].emap[5]].x +
2107 w * nodes[elements[i].emap[6]].x -
2108 v * nodes[elements[i].emap[7]].x -
2109 w * nodes[elements[i].emap[8]].x)) *
2110 (4 * w * nodes[elements[i].emap[9]].y - nodes[elements[i].emap[2]].y +
2111 4 * v * nodes[elements[i].emap[2]].y +
2112 4 * t * nodes[elements[i].emap[5]].y +
2113 4 * u * nodes[elements[i].emap[7]].y) +
2114 ((-1 + 4 * t) * nodes[elements[i].emap[0]].x -
2115 4 * w * nodes[elements[i].emap[9]].x + nodes[elements[i].emap[2]].x -
2116 4 * v * nodes[elements[i].emap[2]].x +
2117 4 * u * nodes[elements[i].emap[4]].x -
2118 4 * t * nodes[elements[i].emap[5]].x +
2119 4 * v * nodes[elements[i].emap[5]].x +
2120 4 * w * nodes[elements[i].emap[6]].x -
2121 4 * u * nodes[elements[i].emap[7]].x) *
2122 ((-1 + 4 * u) * nodes[elements[i].emap[1]].y +
2123 4 * (t * nodes[elements[i].emap[4]].y +
2124 v * nodes[elements[i].emap[7]].y +
2125 w * nodes[elements[i].emap[8]].y));
2126}

Referenced by Coordinates13().

◆ Jacobian3()

void Garfield::ComponentFieldMap::Jacobian3 ( int  i,
double  u,
double  v,
double  w,
double &  det,
double  jac[4][4] 
)
protected

Definition at line 558 of file ComponentFieldMap.cc.

559 {
560
561 // Initial values
562 det = 0;
563 jac[0][0] = 0;
564 jac[0][1] = 0;
565 jac[1][0] = 0;
566 jac[1][1] = 0;
567
568 // Be sure that the element is within range
569 if (i < 0 || i >= nElements) {
570 std::cerr << m_className << "::Jacobian3:\n";
571 std::cerr << " Element " << i << " out of range.\n";
572 return;
573 }
574
575 // Determinant of the quadratic triangular Jacobian
576 det =
577 -(((-1 + 4 * v) * nodes[elements[i].emap[1]].x +
578 nodes[elements[i].emap[2]].x - 4 * w * nodes[elements[i].emap[2]].x +
579 4 * u * nodes[elements[i].emap[3]].x -
580 4 * u * nodes[elements[i].emap[4]].x -
581 4 * v * nodes[elements[i].emap[5]].x +
582 4 * w * nodes[elements[i].emap[5]].x) *
583 (-nodes[elements[i].emap[0]].y + 4 * u * nodes[elements[i].emap[0]].y +
584 4 * v * nodes[elements[i].emap[3]].y +
585 4 * w * nodes[elements[i].emap[4]].y)) -
586 ((-1 + 4 * u) * nodes[elements[i].emap[0]].x +
587 nodes[elements[i].emap[1]].x - 4 * v * nodes[elements[i].emap[1]].x -
588 4 * u * nodes[elements[i].emap[3]].x +
589 4 * v * nodes[elements[i].emap[3]].x +
590 4 * w * nodes[elements[i].emap[4]].x -
591 4 * w * nodes[elements[i].emap[5]].x) *
592 (-nodes[elements[i].emap[2]].y +
593 4 * w * nodes[elements[i].emap[2]].y +
594 4 * u * nodes[elements[i].emap[4]].y +
595 4 * v * nodes[elements[i].emap[5]].y) +
596 ((-1 + 4 * u) * nodes[elements[i].emap[0]].x +
597 nodes[elements[i].emap[2]].x - 4 * w * nodes[elements[i].emap[2]].x +
598 4 * v * nodes[elements[i].emap[3]].x -
599 4 * u * nodes[elements[i].emap[4]].x +
600 4 * w * nodes[elements[i].emap[4]].x -
601 4 * v * nodes[elements[i].emap[5]].x) *
602 (-nodes[elements[i].emap[1]].y +
603 4 * v * nodes[elements[i].emap[1]].y +
604 4 * u * nodes[elements[i].emap[3]].y +
605 4 * w * nodes[elements[i].emap[5]].y);
606
607 // Terms of the quadratic triangular Jacobian
608 jac[0][0] =
609 (-nodes[elements[i].emap[1]].x + 4 * v * nodes[elements[i].emap[1]].x +
610 4 * u * nodes[elements[i].emap[3]].x +
611 4 * w * nodes[elements[i].emap[5]].x) *
612 (-nodes[elements[i].emap[2]].y +
613 4 * w * nodes[elements[i].emap[2]].y +
614 4 * u * nodes[elements[i].emap[4]].y +
615 4 * v * nodes[elements[i].emap[5]].y) -
616 (-nodes[elements[i].emap[2]].x + 4 * w * nodes[elements[i].emap[2]].x +
617 4 * u * nodes[elements[i].emap[4]].x +
618 4 * v * nodes[elements[i].emap[5]].x) *
619 (-nodes[elements[i].emap[1]].y +
620 4 * v * nodes[elements[i].emap[1]].y +
621 4 * u * nodes[elements[i].emap[3]].y +
622 4 * w * nodes[elements[i].emap[5]].y);
623 jac[0][1] = (-1 + 4 * v) * nodes[elements[i].emap[1]].y +
624 nodes[elements[i].emap[2]].y -
625 4 * w * nodes[elements[i].emap[2]].y +
626 4 * u * nodes[elements[i].emap[3]].y -
627 4 * u * nodes[elements[i].emap[4]].y -
628 4 * v * nodes[elements[i].emap[5]].y +
629 4 * w * nodes[elements[i].emap[5]].y;
630 jac[0][2] = nodes[elements[i].emap[1]].x -
631 4 * v * nodes[elements[i].emap[1]].x +
632 (-1 + 4 * w) * nodes[elements[i].emap[2]].x -
633 4 * u * nodes[elements[i].emap[3]].x +
634 4 * u * nodes[elements[i].emap[4]].x +
635 4 * v * nodes[elements[i].emap[5]].x -
636 4 * w * nodes[elements[i].emap[5]].x;
637 jac[1][0] =
638 (-nodes[elements[i].emap[2]].x + 4 * w * nodes[elements[i].emap[2]].x +
639 4 * u * nodes[elements[i].emap[4]].x +
640 4 * v * nodes[elements[i].emap[5]].x) *
641 (-nodes[elements[i].emap[0]].y +
642 4 * u * nodes[elements[i].emap[0]].y +
643 4 * v * nodes[elements[i].emap[3]].y +
644 4 * w * nodes[elements[i].emap[4]].y) -
645 (-nodes[elements[i].emap[0]].x + 4 * u * nodes[elements[i].emap[0]].x +
646 4 * v * nodes[elements[i].emap[3]].x +
647 4 * w * nodes[elements[i].emap[4]].x) *
648 (-nodes[elements[i].emap[2]].y +
649 4 * w * nodes[elements[i].emap[2]].y +
650 4 * u * nodes[elements[i].emap[4]].y +
651 4 * v * nodes[elements[i].emap[5]].y);
652 jac[1][1] =
653 nodes[elements[i].emap[0]].y - 4 * u * nodes[elements[i].emap[0]].y -
654 nodes[elements[i].emap[2]].y + 4 * w * nodes[elements[i].emap[2]].y -
655 4 * v * nodes[elements[i].emap[3]].y +
656 4 * u * nodes[elements[i].emap[4]].y -
657 4 * w * nodes[elements[i].emap[4]].y +
658 4 * v * nodes[elements[i].emap[5]].y;
659 jac[1][2] = (-1 + 4 * u) * nodes[elements[i].emap[0]].x +
660 nodes[elements[i].emap[2]].x -
661 4 * w * nodes[elements[i].emap[2]].x +
662 4 * v * nodes[elements[i].emap[3]].x -
663 4 * u * nodes[elements[i].emap[4]].x +
664 4 * w * nodes[elements[i].emap[4]].x -
665 4 * v * nodes[elements[i].emap[5]].x;
666 jac[2][0] =
667 -((-nodes[elements[i].emap[1]].x + 4 * v * nodes[elements[i].emap[1]].x +
668 4 * u * nodes[elements[i].emap[3]].x +
669 4 * w * nodes[elements[i].emap[5]].x) *
670 (-nodes[elements[i].emap[0]].y + 4 * u * nodes[elements[i].emap[0]].y +
671 4 * v * nodes[elements[i].emap[3]].y +
672 4 * w * nodes[elements[i].emap[4]].y)) +
673 (-nodes[elements[i].emap[0]].x + 4 * u * nodes[elements[i].emap[0]].x +
674 4 * v * nodes[elements[i].emap[3]].x +
675 4 * w * nodes[elements[i].emap[4]].x) *
676 (-nodes[elements[i].emap[1]].y +
677 4 * v * nodes[elements[i].emap[1]].y +
678 4 * u * nodes[elements[i].emap[3]].y +
679 4 * w * nodes[elements[i].emap[5]].y);
680 jac[2][1] = (-1 + 4 * u) * nodes[elements[i].emap[0]].y +
681 nodes[elements[i].emap[1]].y -
682 4 * v * nodes[elements[i].emap[1]].y -
683 4 * u * nodes[elements[i].emap[3]].y +
684 4 * v * nodes[elements[i].emap[3]].y +
685 4 * w * nodes[elements[i].emap[4]].y -
686 4 * w * nodes[elements[i].emap[5]].y;
687 jac[2][2] =
688 nodes[elements[i].emap[0]].x - 4 * u * nodes[elements[i].emap[0]].x -
689 nodes[elements[i].emap[1]].x + 4 * v * nodes[elements[i].emap[1]].x +
690 4 * u * nodes[elements[i].emap[3]].x -
691 4 * v * nodes[elements[i].emap[3]].x -
692 4 * w * nodes[elements[i].emap[4]].x +
693 4 * w * nodes[elements[i].emap[5]].x;
694}

Referenced by Coordinates3().

◆ Jacobian5()

void Garfield::ComponentFieldMap::Jacobian5 ( int  i,
double  u,
double  v,
double &  det,
double  jac[4][4] 
)
protected

Definition at line 696 of file ComponentFieldMap.cc.

697 {
698
699 // Initial values
700 det = 0;
701 jac[0][0] = 0;
702 jac[0][1] = 0;
703 jac[1][0] = 0;
704 jac[1][1] = 0;
705
706 // Be sure that the element is within range
707 if (i < 0 || i >= nElements) {
708 std::cerr << m_className << "::Jacobian5:\n";
709 std::cerr << " Element " << i << " out of range.\n";
710 return;
711 }
712
713 // Determinant of the quadrilateral serendipity Jacobian
714 det =
715 (-2 * u * u * u *
716 ((nodes[elements[i].emap[2]].x + nodes[elements[i].emap[3]].x -
717 2 * nodes[elements[i].emap[6]].x) *
718 (nodes[elements[i].emap[0]].y + nodes[elements[i].emap[1]].y -
719 2 * nodes[elements[i].emap[4]].y) -
720 (nodes[elements[i].emap[0]].x + nodes[elements[i].emap[1]].x -
721 2 * nodes[elements[i].emap[4]].x) *
722 (nodes[elements[i].emap[2]].y + nodes[elements[i].emap[3]].y -
723 2 * nodes[elements[i].emap[6]].y)) +
724 2 * v * v * v *
725 (-((nodes[elements[i].emap[0]].x + nodes[elements[i].emap[3]].x -
726 2 * nodes[elements[i].emap[7]].x) *
727 (nodes[elements[i].emap[1]].y + nodes[elements[i].emap[2]].y -
728 2 * nodes[elements[i].emap[5]].y)) +
729 (nodes[elements[i].emap[1]].x + nodes[elements[i].emap[2]].x -
730 2 * nodes[elements[i].emap[5]].x) *
731 (nodes[elements[i].emap[0]].y + nodes[elements[i].emap[3]].y -
732 2 * nodes[elements[i].emap[7]].y)) +
733 2 * (-((nodes[elements[i].emap[5]].x - nodes[elements[i].emap[7]].x) *
734 (nodes[elements[i].emap[4]].y - nodes[elements[i].emap[6]].y)) +
735 (nodes[elements[i].emap[4]].x - nodes[elements[i].emap[6]].x) *
736 (nodes[elements[i].emap[5]].y - nodes[elements[i].emap[7]].y)) +
737 v * (-(nodes[elements[i].emap[6]].x * nodes[elements[i].emap[0]].y) -
738 2 * nodes[elements[i].emap[7]].x * nodes[elements[i].emap[0]].y +
739 nodes[elements[i].emap[6]].x * nodes[elements[i].emap[1]].y -
740 2 * nodes[elements[i].emap[7]].x * nodes[elements[i].emap[1]].y -
741 nodes[elements[i].emap[6]].x * nodes[elements[i].emap[2]].y -
742 2 * nodes[elements[i].emap[7]].x * nodes[elements[i].emap[2]].y +
743 nodes[elements[i].emap[4]].x *
744 (nodes[elements[i].emap[0]].y - nodes[elements[i].emap[1]].y +
745 nodes[elements[i].emap[2]].y - nodes[elements[i].emap[3]].y) +
746 nodes[elements[i].emap[6]].x * nodes[elements[i].emap[3]].y -
747 2 * nodes[elements[i].emap[7]].x * nodes[elements[i].emap[3]].y -
748 nodes[elements[i].emap[0]].x * nodes[elements[i].emap[4]].y +
749 nodes[elements[i].emap[1]].x * nodes[elements[i].emap[4]].y -
750 nodes[elements[i].emap[2]].x * nodes[elements[i].emap[4]].y +
751 nodes[elements[i].emap[3]].x * nodes[elements[i].emap[4]].y -
752 2 * nodes[elements[i].emap[0]].x * nodes[elements[i].emap[5]].y -
753 2 * nodes[elements[i].emap[1]].x * nodes[elements[i].emap[5]].y -
754 2 * nodes[elements[i].emap[2]].x * nodes[elements[i].emap[5]].y -
755 2 * nodes[elements[i].emap[3]].x * nodes[elements[i].emap[5]].y +
756 8 * nodes[elements[i].emap[7]].x * nodes[elements[i].emap[5]].y +
757 nodes[elements[i].emap[0]].x * nodes[elements[i].emap[6]].y -
758 nodes[elements[i].emap[1]].x * nodes[elements[i].emap[6]].y +
759 nodes[elements[i].emap[2]].x * nodes[elements[i].emap[6]].y -
760 nodes[elements[i].emap[3]].x * nodes[elements[i].emap[6]].y +
761 2 * nodes[elements[i].emap[5]].x *
762 (nodes[elements[i].emap[0]].y + nodes[elements[i].emap[1]].y +
763 nodes[elements[i].emap[2]].y + nodes[elements[i].emap[3]].y -
764 4 * nodes[elements[i].emap[7]].y) +
765 2 * (nodes[elements[i].emap[0]].x + nodes[elements[i].emap[1]].x +
766 nodes[elements[i].emap[2]].x + nodes[elements[i].emap[3]].x) *
767 nodes[elements[i].emap[7]].y) +
768 v * v *
769 (-(nodes[elements[i].emap[4]].x * nodes[elements[i].emap[0]].y) +
770 2 * nodes[elements[i].emap[5]].x * nodes[elements[i].emap[0]].y +
771 nodes[elements[i].emap[6]].x * nodes[elements[i].emap[0]].y +
772 2 * nodes[elements[i].emap[7]].x * nodes[elements[i].emap[0]].y +
773 nodes[elements[i].emap[4]].x * nodes[elements[i].emap[1]].y -
774 2 * nodes[elements[i].emap[5]].x * nodes[elements[i].emap[1]].y -
775 nodes[elements[i].emap[6]].x * nodes[elements[i].emap[1]].y -
776 2 * nodes[elements[i].emap[7]].x * nodes[elements[i].emap[1]].y +
777 nodes[elements[i].emap[4]].x * nodes[elements[i].emap[2]].y +
778 2 * nodes[elements[i].emap[5]].x * nodes[elements[i].emap[2]].y -
779 nodes[elements[i].emap[6]].x * nodes[elements[i].emap[2]].y +
780 2 * nodes[elements[i].emap[7]].x * nodes[elements[i].emap[2]].y -
781 nodes[elements[i].emap[4]].x * nodes[elements[i].emap[3]].y -
782 2 * nodes[elements[i].emap[5]].x * nodes[elements[i].emap[3]].y +
783 nodes[elements[i].emap[6]].x * nodes[elements[i].emap[3]].y -
784 2 * nodes[elements[i].emap[7]].x * nodes[elements[i].emap[3]].y +
785 2 * nodes[elements[i].emap[2]].x *
786 (nodes[elements[i].emap[1]].y + nodes[elements[i].emap[3]].y) -
787 nodes[elements[i].emap[2]].x * nodes[elements[i].emap[4]].y +
788 2 * nodes[elements[i].emap[5]].x * nodes[elements[i].emap[4]].y -
789 2 * nodes[elements[i].emap[7]].x * nodes[elements[i].emap[4]].y -
790 2 * nodes[elements[i].emap[2]].x * nodes[elements[i].emap[5]].y -
791 2 * nodes[elements[i].emap[4]].x * nodes[elements[i].emap[5]].y +
792 2 * nodes[elements[i].emap[6]].x * nodes[elements[i].emap[5]].y +
793 nodes[elements[i].emap[2]].x * nodes[elements[i].emap[6]].y -
794 2 * nodes[elements[i].emap[5]].x * nodes[elements[i].emap[6]].y +
795 2 * nodes[elements[i].emap[7]].x * nodes[elements[i].emap[6]].y +
796 nodes[elements[i].emap[0]].x * (2 * nodes[elements[i].emap[1]].y +
797 2 * nodes[elements[i].emap[3]].y +
798 nodes[elements[i].emap[4]].y -
799 2 * nodes[elements[i].emap[5]].y -
800 nodes[elements[i].emap[6]].y -
801 2 * nodes[elements[i].emap[7]].y) -
802 2 * (nodes[elements[i].emap[2]].x - nodes[elements[i].emap[4]].x +
803 nodes[elements[i].emap[6]].x) *
804 nodes[elements[i].emap[7]].y +
805 nodes[elements[i].emap[3]].x * (-2 * nodes[elements[i].emap[0]].y -
806 2 * nodes[elements[i].emap[2]].y +
807 nodes[elements[i].emap[4]].y +
808 2 * nodes[elements[i].emap[5]].y -
809 nodes[elements[i].emap[6]].y +
810 2 * nodes[elements[i].emap[7]].y) +
811 nodes[elements[i].emap[1]].x * (-2 * nodes[elements[i].emap[0]].y -
812 2 * nodes[elements[i].emap[2]].y -
813 nodes[elements[i].emap[4]].y +
814 2 * nodes[elements[i].emap[5]].y +
815 nodes[elements[i].emap[6]].y +
816 2 * nodes[elements[i].emap[7]].y)) +
817 u * (nodes[elements[i].emap[5]].x * nodes[elements[i].emap[0]].y -
818 2 * nodes[elements[i].emap[6]].x * nodes[elements[i].emap[0]].y -
819 nodes[elements[i].emap[7]].x * nodes[elements[i].emap[0]].y -
820 nodes[elements[i].emap[5]].x * nodes[elements[i].emap[1]].y -
821 2 * nodes[elements[i].emap[6]].x * nodes[elements[i].emap[1]].y +
822 nodes[elements[i].emap[7]].x * nodes[elements[i].emap[1]].y +
823 nodes[elements[i].emap[5]].x * nodes[elements[i].emap[2]].y -
824 2 * nodes[elements[i].emap[6]].x * nodes[elements[i].emap[2]].y -
825 nodes[elements[i].emap[7]].x * nodes[elements[i].emap[2]].y -
826 nodes[elements[i].emap[5]].x * nodes[elements[i].emap[3]].y -
827 2 * nodes[elements[i].emap[6]].x * nodes[elements[i].emap[3]].y +
828 nodes[elements[i].emap[7]].x * nodes[elements[i].emap[3]].y -
829 2 * nodes[elements[i].emap[1]].x * nodes[elements[i].emap[4]].y -
830 2 * nodes[elements[i].emap[2]].x * nodes[elements[i].emap[4]].y -
831 2 * nodes[elements[i].emap[3]].x * nodes[elements[i].emap[4]].y +
832 8 * nodes[elements[i].emap[6]].x * nodes[elements[i].emap[4]].y +
833 nodes[elements[i].emap[1]].x * nodes[elements[i].emap[5]].y -
834 nodes[elements[i].emap[2]].x * nodes[elements[i].emap[5]].y +
835 nodes[elements[i].emap[3]].x * nodes[elements[i].emap[5]].y +
836 2 * nodes[elements[i].emap[4]].x *
837 (nodes[elements[i].emap[0]].y + nodes[elements[i].emap[1]].y +
838 nodes[elements[i].emap[2]].y + nodes[elements[i].emap[3]].y -
839 4 * nodes[elements[i].emap[6]].y) +
840 2 * nodes[elements[i].emap[1]].x * nodes[elements[i].emap[6]].y +
841 2 * nodes[elements[i].emap[2]].x * nodes[elements[i].emap[6]].y +
842 2 * nodes[elements[i].emap[3]].x * nodes[elements[i].emap[6]].y -
843 (nodes[elements[i].emap[1]].x - nodes[elements[i].emap[2]].x +
844 nodes[elements[i].emap[3]].x) *
845 nodes[elements[i].emap[7]].y +
846 nodes[elements[i].emap[0]].x * (-2 * nodes[elements[i].emap[4]].y -
847 nodes[elements[i].emap[5]].y +
848 2 * nodes[elements[i].emap[6]].y +
849 nodes[elements[i].emap[7]].y) +
850 v * v * (4 * nodes[elements[i].emap[4]].x *
851 nodes[elements[i].emap[0]].y -
852 3 * nodes[elements[i].emap[5]].x *
853 nodes[elements[i].emap[0]].y -
854 4 * nodes[elements[i].emap[6]].x *
855 nodes[elements[i].emap[0]].y -
856 5 * nodes[elements[i].emap[7]].x *
857 nodes[elements[i].emap[0]].y +
858 4 * nodes[elements[i].emap[4]].x *
859 nodes[elements[i].emap[1]].y -
860 5 * nodes[elements[i].emap[5]].x *
861 nodes[elements[i].emap[1]].y -
862 4 * nodes[elements[i].emap[6]].x *
863 nodes[elements[i].emap[1]].y -
864 3 * nodes[elements[i].emap[7]].x *
865 nodes[elements[i].emap[1]].y +
866 4 * nodes[elements[i].emap[4]].x *
867 nodes[elements[i].emap[2]].y +
868 5 * nodes[elements[i].emap[5]].x *
869 nodes[elements[i].emap[2]].y -
870 4 * nodes[elements[i].emap[6]].x *
871 nodes[elements[i].emap[2]].y +
872 3 * nodes[elements[i].emap[7]].x *
873 nodes[elements[i].emap[2]].y +
874 4 * nodes[elements[i].emap[4]].x *
875 nodes[elements[i].emap[3]].y +
876 3 * nodes[elements[i].emap[5]].x *
877 nodes[elements[i].emap[3]].y -
878 4 * nodes[elements[i].emap[6]].x *
879 nodes[elements[i].emap[3]].y +
880 5 * nodes[elements[i].emap[7]].x *
881 nodes[elements[i].emap[3]].y +
882 8 * nodes[elements[i].emap[5]].x *
883 nodes[elements[i].emap[4]].y +
884 8 * nodes[elements[i].emap[7]].x *
885 nodes[elements[i].emap[4]].y -
886 8 * nodes[elements[i].emap[4]].x *
887 nodes[elements[i].emap[5]].y +
888 8 * nodes[elements[i].emap[6]].x *
889 nodes[elements[i].emap[5]].y -
890 8 * nodes[elements[i].emap[5]].x *
891 nodes[elements[i].emap[6]].y -
892 8 * nodes[elements[i].emap[7]].x *
893 nodes[elements[i].emap[6]].y +
894 nodes[elements[i].emap[3]].x *
895 (5 * nodes[elements[i].emap[0]].y +
896 3 * nodes[elements[i].emap[1]].y -
897 4 * nodes[elements[i].emap[4]].y -
898 3 * nodes[elements[i].emap[5]].y +
899 4 * nodes[elements[i].emap[6]].y -
900 5 * nodes[elements[i].emap[7]].y) +
901 nodes[elements[i].emap[2]].x *
902 (3 * nodes[elements[i].emap[0]].y +
903 5 * nodes[elements[i].emap[1]].y -
904 4 * nodes[elements[i].emap[4]].y -
905 5 * nodes[elements[i].emap[5]].y +
906 4 * nodes[elements[i].emap[6]].y -
907 3 * nodes[elements[i].emap[7]].y) -
908 8 * nodes[elements[i].emap[4]].x *
909 nodes[elements[i].emap[7]].y +
910 8 * nodes[elements[i].emap[6]].x *
911 nodes[elements[i].emap[7]].y +
912 nodes[elements[i].emap[1]].x *
913 (-5 * nodes[elements[i].emap[2]].y -
914 3 * nodes[elements[i].emap[3]].y -
915 4 * nodes[elements[i].emap[4]].y +
916 5 * nodes[elements[i].emap[5]].y +
917 4 * nodes[elements[i].emap[6]].y +
918 3 * nodes[elements[i].emap[7]].y) +
919 nodes[elements[i].emap[0]].x *
920 (-3 * nodes[elements[i].emap[2]].y -
921 5 * nodes[elements[i].emap[3]].y -
922 4 * nodes[elements[i].emap[4]].y +
923 3 * nodes[elements[i].emap[5]].y +
924 4 * nodes[elements[i].emap[6]].y +
925 5 * nodes[elements[i].emap[7]].y)) -
926 2 * v *
927 (nodes[elements[i].emap[6]].x * nodes[elements[i].emap[0]].y -
928 3 * nodes[elements[i].emap[7]].x *
929 nodes[elements[i].emap[0]].y +
930 nodes[elements[i].emap[6]].x * nodes[elements[i].emap[1]].y -
931 nodes[elements[i].emap[7]].x * nodes[elements[i].emap[1]].y +
932 3 * nodes[elements[i].emap[6]].x *
933 nodes[elements[i].emap[2]].y -
934 nodes[elements[i].emap[7]].x * nodes[elements[i].emap[2]].y +
935 3 * nodes[elements[i].emap[6]].x *
936 nodes[elements[i].emap[3]].y -
937 3 * nodes[elements[i].emap[7]].x *
938 nodes[elements[i].emap[3]].y -
939 3 * nodes[elements[i].emap[0]].x *
940 nodes[elements[i].emap[4]].y -
941 3 * nodes[elements[i].emap[1]].x *
942 nodes[elements[i].emap[4]].y -
943 nodes[elements[i].emap[2]].x * nodes[elements[i].emap[4]].y -
944 nodes[elements[i].emap[3]].x * nodes[elements[i].emap[4]].y +
945 4 * nodes[elements[i].emap[7]].x *
946 nodes[elements[i].emap[4]].y +
947 nodes[elements[i].emap[0]].x * nodes[elements[i].emap[5]].y +
948 3 * nodes[elements[i].emap[1]].x *
949 nodes[elements[i].emap[5]].y +
950 3 * nodes[elements[i].emap[2]].x *
951 nodes[elements[i].emap[5]].y +
952 nodes[elements[i].emap[3]].x * nodes[elements[i].emap[5]].y -
953 4 * nodes[elements[i].emap[6]].x *
954 nodes[elements[i].emap[5]].y -
955 nodes[elements[i].emap[0]].x * nodes[elements[i].emap[6]].y -
956 nodes[elements[i].emap[1]].x * nodes[elements[i].emap[6]].y -
957 3 * nodes[elements[i].emap[2]].x *
958 nodes[elements[i].emap[6]].y -
959 3 * nodes[elements[i].emap[3]].x *
960 nodes[elements[i].emap[6]].y +
961 4 * nodes[elements[i].emap[7]].x *
962 nodes[elements[i].emap[6]].y -
963 nodes[elements[i].emap[5]].x *
964 (nodes[elements[i].emap[0]].y +
965 3 * nodes[elements[i].emap[1]].y +
966 3 * nodes[elements[i].emap[2]].y +
967 nodes[elements[i].emap[3]].y -
968 4 * (nodes[elements[i].emap[4]].y +
969 nodes[elements[i].emap[6]].y)) +
970 (3 * nodes[elements[i].emap[0]].x +
971 nodes[elements[i].emap[1]].x + nodes[elements[i].emap[2]].x +
972 3 * nodes[elements[i].emap[3]].x -
973 4 * nodes[elements[i].emap[6]].x) *
974 nodes[elements[i].emap[7]].y +
975 nodes[elements[i].emap[4]].x *
976 (3 * nodes[elements[i].emap[0]].y +
977 3 * nodes[elements[i].emap[1]].y +
978 nodes[elements[i].emap[2]].y +
979 nodes[elements[i].emap[3]].y -
980 4 * (nodes[elements[i].emap[5]].y +
981 nodes[elements[i].emap[7]].y)))) +
982 u * u *
983 (2 * nodes[elements[i].emap[3]].x * nodes[elements[i].emap[0]].y -
984 2 * nodes[elements[i].emap[4]].x * nodes[elements[i].emap[0]].y -
985 nodes[elements[i].emap[5]].x * nodes[elements[i].emap[0]].y -
986 2 * nodes[elements[i].emap[6]].x * nodes[elements[i].emap[0]].y +
987 nodes[elements[i].emap[7]].x * nodes[elements[i].emap[0]].y -
988 2 * nodes[elements[i].emap[0]].x * nodes[elements[i].emap[1]].y +
989 2 * nodes[elements[i].emap[4]].x * nodes[elements[i].emap[1]].y -
990 nodes[elements[i].emap[5]].x * nodes[elements[i].emap[1]].y +
991 2 * nodes[elements[i].emap[6]].x * nodes[elements[i].emap[1]].y +
992 nodes[elements[i].emap[7]].x * nodes[elements[i].emap[1]].y +
993 2 * nodes[elements[i].emap[3]].x * nodes[elements[i].emap[2]].y -
994 2 * nodes[elements[i].emap[4]].x * nodes[elements[i].emap[2]].y +
995 nodes[elements[i].emap[5]].x * nodes[elements[i].emap[2]].y -
996 2 * nodes[elements[i].emap[6]].x * nodes[elements[i].emap[2]].y -
997 nodes[elements[i].emap[7]].x * nodes[elements[i].emap[2]].y +
998 2 * nodes[elements[i].emap[4]].x * nodes[elements[i].emap[3]].y +
999 nodes[elements[i].emap[5]].x * nodes[elements[i].emap[3]].y +
1000 2 * nodes[elements[i].emap[6]].x * nodes[elements[i].emap[3]].y -
1001 nodes[elements[i].emap[7]].x * nodes[elements[i].emap[3]].y -
1002 2 * nodes[elements[i].emap[3]].x * nodes[elements[i].emap[4]].y +
1003 2 * nodes[elements[i].emap[5]].x * nodes[elements[i].emap[4]].y -
1004 2 * nodes[elements[i].emap[7]].x * nodes[elements[i].emap[4]].y -
1005 nodes[elements[i].emap[3]].x * nodes[elements[i].emap[5]].y -
1006 2 * nodes[elements[i].emap[4]].x * nodes[elements[i].emap[5]].y +
1007 2 * nodes[elements[i].emap[6]].x * nodes[elements[i].emap[5]].y -
1008 2 * nodes[elements[i].emap[3]].x * nodes[elements[i].emap[6]].y -
1009 2 * nodes[elements[i].emap[5]].x * nodes[elements[i].emap[6]].y +
1010 2 * nodes[elements[i].emap[7]].x * nodes[elements[i].emap[6]].y +
1011 nodes[elements[i].emap[0]].x * (-2 * nodes[elements[i].emap[3]].y +
1012 2 * nodes[elements[i].emap[4]].y +
1013 nodes[elements[i].emap[5]].y +
1014 2 * nodes[elements[i].emap[6]].y -
1015 nodes[elements[i].emap[7]].y) +
1016 (nodes[elements[i].emap[3]].x + 2 * nodes[elements[i].emap[4]].x -
1017 2 * nodes[elements[i].emap[6]].x) *
1018 nodes[elements[i].emap[7]].y +
1019 nodes[elements[i].emap[2]].x * (-2 * nodes[elements[i].emap[1]].y -
1020 2 * nodes[elements[i].emap[3]].y +
1021 2 * nodes[elements[i].emap[4]].y -
1022 nodes[elements[i].emap[5]].y +
1023 2 * nodes[elements[i].emap[6]].y +
1024 nodes[elements[i].emap[7]].y) -
1025 3 * v * v *
1026 (nodes[elements[i].emap[5]].x * nodes[elements[i].emap[0]].y -
1027 nodes[elements[i].emap[6]].x * nodes[elements[i].emap[0]].y -
1028 nodes[elements[i].emap[7]].x * nodes[elements[i].emap[0]].y +
1029 nodes[elements[i].emap[5]].x * nodes[elements[i].emap[1]].y +
1030 nodes[elements[i].emap[6]].x * nodes[elements[i].emap[1]].y -
1031 nodes[elements[i].emap[7]].x * nodes[elements[i].emap[1]].y -
1032 nodes[elements[i].emap[5]].x * nodes[elements[i].emap[2]].y +
1033 nodes[elements[i].emap[6]].x * nodes[elements[i].emap[2]].y +
1034 nodes[elements[i].emap[7]].x * nodes[elements[i].emap[2]].y -
1035 nodes[elements[i].emap[5]].x * nodes[elements[i].emap[3]].y -
1036 nodes[elements[i].emap[6]].x * nodes[elements[i].emap[3]].y +
1037 nodes[elements[i].emap[7]].x * nodes[elements[i].emap[3]].y -
1038 2 * nodes[elements[i].emap[5]].x *
1039 nodes[elements[i].emap[4]].y +
1040 2 * nodes[elements[i].emap[7]].x *
1041 nodes[elements[i].emap[4]].y -
1042 2 * nodes[elements[i].emap[6]].x *
1043 nodes[elements[i].emap[5]].y +
1044 2 * nodes[elements[i].emap[5]].x *
1045 nodes[elements[i].emap[6]].y -
1046 2 * nodes[elements[i].emap[7]].x *
1047 nodes[elements[i].emap[6]].y +
1048 nodes[elements[i].emap[4]].x *
1049 (nodes[elements[i].emap[0]].y -
1050 nodes[elements[i].emap[1]].y -
1051 nodes[elements[i].emap[2]].y +
1052 nodes[elements[i].emap[3]].y +
1053 2 * nodes[elements[i].emap[5]].y -
1054 2 * nodes[elements[i].emap[7]].y) +
1055 nodes[elements[i].emap[3]].x * (nodes[elements[i].emap[0]].y -
1056 nodes[elements[i].emap[2]].y -
1057 nodes[elements[i].emap[4]].y +
1058 nodes[elements[i].emap[5]].y +
1059 nodes[elements[i].emap[6]].y -
1060 nodes[elements[i].emap[7]].y) +
1061 2 * nodes[elements[i].emap[6]].x *
1062 nodes[elements[i].emap[7]].y +
1063 (nodes[elements[i].emap[0]].x - nodes[elements[i].emap[2]].x) *
1064 (nodes[elements[i].emap[1]].y -
1065 nodes[elements[i].emap[3]].y -
1066 nodes[elements[i].emap[4]].y -
1067 nodes[elements[i].emap[5]].y +
1068 nodes[elements[i].emap[6]].y +
1069 nodes[elements[i].emap[7]].y)) +
1070 v * (4 * nodes[elements[i].emap[5]].x *
1071 nodes[elements[i].emap[0]].y +
1072 3 * nodes[elements[i].emap[6]].x *
1073 nodes[elements[i].emap[0]].y -
1074 4 * nodes[elements[i].emap[7]].x *
1075 nodes[elements[i].emap[0]].y +
1076 4 * nodes[elements[i].emap[5]].x *
1077 nodes[elements[i].emap[1]].y -
1078 3 * nodes[elements[i].emap[6]].x *
1079 nodes[elements[i].emap[1]].y -
1080 4 * nodes[elements[i].emap[7]].x *
1081 nodes[elements[i].emap[1]].y +
1082 4 * nodes[elements[i].emap[5]].x *
1083 nodes[elements[i].emap[2]].y -
1084 5 * nodes[elements[i].emap[6]].x *
1085 nodes[elements[i].emap[2]].y -
1086 4 * nodes[elements[i].emap[7]].x *
1087 nodes[elements[i].emap[2]].y +
1088 4 * nodes[elements[i].emap[5]].x *
1089 nodes[elements[i].emap[3]].y +
1090 5 * nodes[elements[i].emap[6]].x *
1091 nodes[elements[i].emap[3]].y -
1092 4 * nodes[elements[i].emap[7]].x *
1093 nodes[elements[i].emap[3]].y -
1094 8 * nodes[elements[i].emap[5]].x *
1095 nodes[elements[i].emap[4]].y +
1096 8 * nodes[elements[i].emap[7]].x *
1097 nodes[elements[i].emap[4]].y +
1098 8 * nodes[elements[i].emap[6]].x *
1099 nodes[elements[i].emap[5]].y -
1100 8 * nodes[elements[i].emap[5]].x *
1101 nodes[elements[i].emap[6]].y +
1102 8 * nodes[elements[i].emap[7]].x *
1103 nodes[elements[i].emap[6]].y +
1104 nodes[elements[i].emap[4]].x *
1105 (5 * nodes[elements[i].emap[0]].y -
1106 5 * nodes[elements[i].emap[1]].y -
1107 3 * nodes[elements[i].emap[2]].y +
1108 3 * nodes[elements[i].emap[3]].y +
1109 8 * nodes[elements[i].emap[5]].y -
1110 8 * nodes[elements[i].emap[7]].y) -
1111 8 * nodes[elements[i].emap[6]].x *
1112 nodes[elements[i].emap[7]].y +
1113 nodes[elements[i].emap[3]].x *
1114 (3 * nodes[elements[i].emap[1]].y +
1115 5 * nodes[elements[i].emap[2]].y -
1116 3 * nodes[elements[i].emap[4]].y -
1117 4 * nodes[elements[i].emap[5]].y -
1118 5 * nodes[elements[i].emap[6]].y +
1119 4 * nodes[elements[i].emap[7]].y) +
1120 nodes[elements[i].emap[0]].x *
1121 (5 * nodes[elements[i].emap[1]].y +
1122 3 * nodes[elements[i].emap[2]].y -
1123 5 * nodes[elements[i].emap[4]].y -
1124 4 * nodes[elements[i].emap[5]].y -
1125 3 * nodes[elements[i].emap[6]].y +
1126 4 * nodes[elements[i].emap[7]].y) +
1127 nodes[elements[i].emap[2]].x *
1128 (-3 * nodes[elements[i].emap[0]].y -
1129 5 * nodes[elements[i].emap[3]].y +
1130 3 * nodes[elements[i].emap[4]].y -
1131 4 * nodes[elements[i].emap[5]].y +
1132 5 * nodes[elements[i].emap[6]].y +
1133 4 * nodes[elements[i].emap[7]].y)) +
1134 nodes[elements[i].emap[1]].x *
1135 ((-1 + v) * (-2 + 3 * v) * nodes[elements[i].emap[0]].y +
1136 2 * nodes[elements[i].emap[2]].y -
1137 2 * nodes[elements[i].emap[4]].y +
1138 nodes[elements[i].emap[5]].y -
1139 2 * nodes[elements[i].emap[6]].y -
1140 nodes[elements[i].emap[7]].y +
1141 v * (-3 * nodes[elements[i].emap[3]].y +
1142 5 * nodes[elements[i].emap[4]].y -
1143 4 * nodes[elements[i].emap[5]].y +
1144 3 * nodes[elements[i].emap[6]].y +
1145 4 * nodes[elements[i].emap[7]].y -
1146 3 * v * (nodes[elements[i].emap[2]].y +
1147 nodes[elements[i].emap[4]].y -
1148 nodes[elements[i].emap[5]].y -
1149 nodes[elements[i].emap[6]].y +
1150 nodes[elements[i].emap[7]].y))))) /
1151 8;
1152 // Jacobian terms
1153 jac[0][0] =
1154 (u * u * (-nodes[elements[i].emap[0]].y - nodes[elements[i].emap[1]].y +
1155 nodes[elements[i].emap[2]].y + nodes[elements[i].emap[3]].y +
1156 2 * nodes[elements[i].emap[4]].y -
1157 2 * nodes[elements[i].emap[6]].y) +
1158 2 * (-nodes[elements[i].emap[4]].y + nodes[elements[i].emap[6]].y +
1159 v * (nodes[elements[i].emap[0]].y + nodes[elements[i].emap[1]].y +
1160 nodes[elements[i].emap[2]].y + nodes[elements[i].emap[3]].y -
1161 2 * nodes[elements[i].emap[5]].y -
1162 2 * nodes[elements[i].emap[7]].y)) +
1163 u * (nodes[elements[i].emap[0]].y -
1164 2 * v * nodes[elements[i].emap[0]].y -
1165 nodes[elements[i].emap[1]].y +
1166 2 * v * nodes[elements[i].emap[1]].y +
1167 nodes[elements[i].emap[2]].y +
1168 2 * v * nodes[elements[i].emap[2]].y -
1169 nodes[elements[i].emap[3]].y -
1170 2 * v * nodes[elements[i].emap[3]].y -
1171 4 * v * nodes[elements[i].emap[5]].y +
1172 4 * v * nodes[elements[i].emap[7]].y)) /
1173 4;
1174 jac[0][1] =
1175 (u * u * (nodes[elements[i].emap[0]].x + nodes[elements[i].emap[1]].x -
1176 nodes[elements[i].emap[2]].x - nodes[elements[i].emap[3]].x -
1177 2 * nodes[elements[i].emap[4]].x +
1178 2 * nodes[elements[i].emap[6]].x) -
1179 2 * (-nodes[elements[i].emap[4]].x + nodes[elements[i].emap[6]].x +
1180 v * (nodes[elements[i].emap[0]].x + nodes[elements[i].emap[1]].x +
1181 nodes[elements[i].emap[2]].x + nodes[elements[i].emap[3]].x -
1182 2 * nodes[elements[i].emap[5]].x -
1183 2 * nodes[elements[i].emap[7]].x)) +
1184 u * ((-1 + 2 * v) * nodes[elements[i].emap[0]].x +
1185 nodes[elements[i].emap[1]].x -
1186 2 * v * nodes[elements[i].emap[1]].x -
1187 nodes[elements[i].emap[2]].x -
1188 2 * v * nodes[elements[i].emap[2]].x +
1189 nodes[elements[i].emap[3]].x +
1190 2 * v * nodes[elements[i].emap[3]].x +
1191 4 * v * nodes[elements[i].emap[5]].x -
1192 4 * v * nodes[elements[i].emap[7]].x)) /
1193 4;
1194 jac[1][0] =
1195 (v * (-nodes[elements[i].emap[0]].y + nodes[elements[i].emap[1]].y -
1196 nodes[elements[i].emap[2]].y + nodes[elements[i].emap[3]].y) -
1197 2 * nodes[elements[i].emap[5]].y +
1198 2 * u *
1199 ((-1 + v) * nodes[elements[i].emap[0]].y +
1200 (-1 + v) * nodes[elements[i].emap[1]].y -
1201 nodes[elements[i].emap[2]].y - v * nodes[elements[i].emap[2]].y -
1202 nodes[elements[i].emap[3]].y - v * nodes[elements[i].emap[3]].y +
1203 2 * nodes[elements[i].emap[4]].y -
1204 2 * v * nodes[elements[i].emap[4]].y +
1205 2 * nodes[elements[i].emap[6]].y +
1206 2 * v * nodes[elements[i].emap[6]].y) +
1207 v * v * (nodes[elements[i].emap[0]].y - nodes[elements[i].emap[1]].y -
1208 nodes[elements[i].emap[2]].y + nodes[elements[i].emap[3]].y +
1209 2 * nodes[elements[i].emap[5]].y -
1210 2 * nodes[elements[i].emap[7]].y) +
1211 2 * nodes[elements[i].emap[7]].y) /
1212 4;
1213 jac[1][1] =
1214 (v * (nodes[elements[i].emap[0]].x - nodes[elements[i].emap[1]].x +
1215 nodes[elements[i].emap[2]].x - nodes[elements[i].emap[3]].x) +
1216 2 * u *
1217 (nodes[elements[i].emap[0]].x - v * nodes[elements[i].emap[0]].x +
1218 nodes[elements[i].emap[1]].x - v * nodes[elements[i].emap[1]].x +
1219 nodes[elements[i].emap[2]].x + v * nodes[elements[i].emap[2]].x +
1220 nodes[elements[i].emap[3]].x + v * nodes[elements[i].emap[3]].x -
1221 2 * nodes[elements[i].emap[4]].x +
1222 2 * v * nodes[elements[i].emap[4]].x -
1223 2 * nodes[elements[i].emap[6]].x -
1224 2 * v * nodes[elements[i].emap[6]].x) +
1225 2 * (nodes[elements[i].emap[5]].x - nodes[elements[i].emap[7]].x) +
1226 v * v * (-nodes[elements[i].emap[0]].x + nodes[elements[i].emap[1]].x +
1227 nodes[elements[i].emap[2]].x - nodes[elements[i].emap[3]].x -
1228 2 * nodes[elements[i].emap[5]].x +
1229 2 * nodes[elements[i].emap[7]].x)) /
1230 4;
1231}

Referenced by Coordinates5().

◆ JacobianCube()

void Garfield::ComponentFieldMap::JacobianCube ( int  i,
double  t1,
double  t2,
double  t3,
TMatrixD *&  jac,
std::vector< TMatrixD * > &  dN 
)
protected

Definition at line 2128 of file ComponentFieldMap.cc.

2130 {
2131 if (jac == 0) {
2132 std::cerr << m_className << "::JacobianCube:\n";
2133 std::cerr << " Pointer to Jacobian matrix is empty!\n";
2134 return;
2135 }
2136 dN.clear();
2137 // Be sure that the element is within range
2138 if (element < 0 || element >= nElements) {
2139 std::cerr << m_className << "::JacobianCube:\n";
2140 std::cerr << " Element " << element << " out of range.\n";
2141 return;
2142 }
2143 // Here the partial derivatives of the 8 shaping functions are calculated
2144 double N1[3] = {-1 * (1 - t2) * (1 - t3), (1 - t1) * -1 * (1 - t3),
2145 (1 - t1) * (1 - t2) * -1};
2146 double N2[3] = {+1 * (1 - t2) * (1 - t3), (1 + t1) * -1 * (1 - t3),
2147 (1 + t1) * (1 - t2) * -1};
2148 double N3[3] = {+1 * (1 + t2) * (1 - t3), (1 + t1) * +1 * (1 - t3),
2149 (1 + t1) * (1 + t2) * -1};
2150 double N4[3] = {-1 * (1 + t2) * (1 - t3), (1 - t1) * +1 * (1 - t3),
2151 (1 - t1) * (1 + t2) * -1};
2152 double N5[3] = {-1 * (1 - t2) * (1 + t3), (1 - t1) * -1 * (1 + t3),
2153 (1 - t1) * (1 - t2) * +1};
2154 double N6[3] = {+1 * (1 - t2) * (1 + t3), (1 + t1) * -1 * (1 + t3),
2155 (1 + t1) * (1 - t2) * +1};
2156 double N7[3] = {+1 * (1 + t2) * (1 + t3), (1 + t1) * +1 * (1 + t3),
2157 (1 + t1) * (1 + t2) * +1};
2158 double N8[3] = {-1 * (1 + t2) * (1 + t3), (1 - t1) * +1 * (1 + t3),
2159 (1 - t1) * (1 + t2) * +1};
2160 // Partial derivatives are stored in dN
2161 TMatrixD* m_N1 = new TMatrixD(3, 1, N1);
2162 *m_N1 = (1. / 8. * (*m_N1));
2163 dN.push_back(m_N1);
2164 TMatrixD* m_N2 = new TMatrixD(3, 1, N2);
2165 *m_N2 = (1. / 8. * (*m_N2));
2166 dN.push_back(m_N2);
2167 TMatrixD* m_N3 = new TMatrixD(3, 1, N3);
2168 *m_N3 = (1. / 8. * (*m_N3));
2169 dN.push_back(m_N3);
2170 TMatrixD* m_N4 = new TMatrixD(3, 1, N4);
2171 *m_N4 = (1. / 8. * (*m_N4));
2172 dN.push_back(m_N4);
2173 TMatrixD* m_N5 = new TMatrixD(3, 1, N5);
2174 *m_N5 = (1. / 8. * (*m_N5));
2175 dN.push_back(m_N5);
2176 TMatrixD* m_N6 = new TMatrixD(3, 1, N6);
2177 *m_N6 = (1. / 8. * (*m_N6));
2178 dN.push_back(m_N6);
2179 TMatrixD* m_N7 = new TMatrixD(3, 1, N7);
2180 *m_N7 = (1. / 8. * (*m_N7));
2181 dN.push_back(m_N7);
2182 TMatrixD* m_N8 = new TMatrixD(3, 1, N8);
2183 *m_N8 = (1. / 8. * (*m_N8));
2184 dN.push_back(m_N8);
2185 // Calculation of the jacobian using dN
2186 for (int node = 0; node < 8; node++) {
2187 (*jac)(0, 0) +=
2188 nodes[elements[element].emap[node]].x * ((*dN.at(node))(0, 0));
2189 (*jac)(0, 1) +=
2190 nodes[elements[element].emap[node]].y * ((*dN.at(node))(0, 0));
2191 (*jac)(0, 2) +=
2192 nodes[elements[element].emap[node]].z * ((*dN.at(node))(0, 0));
2193 (*jac)(1, 0) +=
2194 nodes[elements[element].emap[node]].x * ((*dN.at(node))(1, 0));
2195 (*jac)(1, 1) +=
2196 nodes[elements[element].emap[node]].y * ((*dN.at(node))(1, 0));
2197 (*jac)(1, 2) +=
2198 nodes[elements[element].emap[node]].z * ((*dN.at(node))(1, 0));
2199 (*jac)(2, 0) +=
2200 nodes[elements[element].emap[node]].x * ((*dN.at(node))(2, 0));
2201 (*jac)(2, 1) +=
2202 nodes[elements[element].emap[node]].y * ((*dN.at(node))(2, 0));
2203 (*jac)(2, 2) +=
2204 nodes[elements[element].emap[node]].z * ((*dN.at(node))(2, 0));
2205 }
2206
2207 // compute determinant
2208 if (debug) {
2209 std::cout << m_className << "::JacobianCube:" << std::endl;
2210 std::cout << " Det.: " << jac->Determinant() << std::endl;
2211 std::cout << " Jacobian matrix.: " << std::endl;
2212 jac->Print("%11.10g");
2213 std::cout << " Hexahedral coordinates (t, u, v) = (" << t1 << "," << t2
2214 << "," << t3 << ")" << std::endl;
2215 std::cout << " Node xyzV" << std::endl;
2216 for (int node = 0; node < 8; node++) {
2217 std::cout << " " << elements[element].emap[node] << " "
2218 << nodes[elements[element].emap[node]].x << " "
2219 << nodes[elements[element].emap[node]].y << " "
2220 << nodes[elements[element].emap[node]].z << " "
2221 << nodes[elements[element].emap[node]].v << std::endl;
2222 }
2223 }
2224}

Referenced by CoordinatesCube().

◆ MapCoordinates()

void Garfield::ComponentFieldMap::MapCoordinates ( double &  xpos,
double &  ypos,
double &  zpos,
bool &  xmirrored,
bool &  ymirrored,
bool &  zmirrored,
double &  rcoordinate,
double &  rotation 
) const
protected

Definition at line 4036 of file ComponentFieldMap.cc.

4039 {
4040
4041 // Initial values
4042 rotation = 0;
4043
4044 // If chamber is periodic, reduce to the cell volume.
4045 xmirrored = false;
4046 double auxr, auxphi;
4047 if (xPeriodic) {
4048 xpos = mapxmin + fmod(xpos - mapxmin, mapxmax - mapxmin);
4049 if (xpos < mapxmin) xpos += mapxmax - mapxmin;
4050 } else if (xMirrorPeriodic) {
4051 double xnew = mapxmin + fmod(xpos - mapxmin, mapxmax - mapxmin);
4052 if (xnew < mapxmin) xnew += mapxmax - mapxmin;
4053 int nx = int(floor(0.5 + (xnew - xpos) / (mapxmax - mapxmin)));
4054 if (nx != 2 * (nx / 2)) {
4055 xnew = mapxmin + mapxmax - xnew;
4056 xmirrored = true;
4057 }
4058 xpos = xnew;
4059 }
4060 if (xAxiallyPeriodic && (zpos != 0 || ypos != 0)) {
4061 auxr = sqrt(zpos * zpos + ypos * ypos);
4062 auxphi = atan2(zpos, ypos);
4063 rotation = (mapxamax - mapxamin) *
4064 floor(0.5 + (auxphi - 0.5 * (mapxamin + mapxamax)) /
4065 (mapxamax - mapxamin));
4066 if (auxphi - rotation < mapxamin)
4067 rotation = rotation - (mapxamax - mapxamin);
4068 if (auxphi - rotation > mapxamax)
4069 rotation = rotation + (mapxamax - mapxamin);
4070 auxphi = auxphi - rotation;
4071 ypos = auxr * cos(auxphi);
4072 zpos = auxr * sin(auxphi);
4073 }
4074
4075 ymirrored = false;
4076 if (yPeriodic) {
4077 ypos = mapymin + fmod(ypos - mapymin, mapymax - mapymin);
4078 if (ypos < mapymin) ypos += mapymax - mapymin;
4079 } else if (yMirrorPeriodic) {
4080 double ynew = mapymin + fmod(ypos - mapymin, mapymax - mapymin);
4081 if (ynew < mapymin) ynew += mapymax - mapymin;
4082 int ny = int(floor(0.5 + (ynew - ypos) / (mapymax - mapymin)));
4083 if (ny != 2 * (ny / 2)) {
4084 ynew = mapymin + mapymax - ynew;
4085 ymirrored = true;
4086 }
4087 ypos = ynew;
4088 }
4089 if (yAxiallyPeriodic && (xpos != 0 || zpos != 0)) {
4090 auxr = sqrt(xpos * xpos + zpos * zpos);
4091 auxphi = atan2(xpos, zpos);
4092 rotation = (mapyamax - mapyamin) *
4093 floor(0.5 + (auxphi - 0.5 * (mapyamin + mapyamax)) /
4094 (mapyamax - mapyamin));
4095 if (auxphi - rotation < mapyamin)
4096 rotation = rotation - (mapyamax - mapyamin);
4097 if (auxphi - rotation > mapyamax)
4098 rotation = rotation + (mapyamax - mapyamin);
4099 auxphi = auxphi - rotation;
4100 zpos = auxr * cos(auxphi);
4101 xpos = auxr * sin(auxphi);
4102 }
4103
4104 zmirrored = false;
4105 if (zPeriodic) {
4106 zpos = mapzmin + fmod(zpos - mapzmin, mapzmax - mapzmin);
4107 if (zpos < mapzmin) zpos += mapzmax - mapzmin;
4108 } else if (zMirrorPeriodic) {
4109 double znew = mapzmin + fmod(zpos - mapzmin, mapzmax - mapzmin);
4110 if (znew < mapzmin) znew += mapzmax - mapzmin;
4111 int nz = int(floor(0.5 + (znew - zpos) / (mapzmax - mapzmin)));
4112 if (nz != 2 * (nz / 2)) {
4113 znew = mapzmin + mapzmax - znew;
4114 zmirrored = true;
4115 }
4116 zpos = znew;
4117 }
4118 if (zAxiallyPeriodic && (ypos != 0 || xpos != 0)) {
4119 auxr = sqrt(ypos * ypos + xpos * xpos);
4120 auxphi = atan2(ypos, xpos);
4121 rotation = (mapzamax - mapzamin) *
4122 floor(0.5 + (auxphi - 0.5 * (mapzamin + mapzamax)) /
4123 (mapzamax - mapzamin));
4124 if (auxphi - rotation < mapzamin)
4125 rotation = rotation - (mapzamax - mapzamin);
4126 if (auxphi - rotation > mapzamax)
4127 rotation = rotation + (mapzamax - mapzamin);
4128 auxphi = auxphi - rotation;
4129 xpos = auxr * cos(auxphi);
4130 ypos = auxr * sin(auxphi);
4131 }
4132
4133 // If we have a rotationally symmetric field map, store coordinates.
4134 rcoordinate = 0;
4135 double zcoordinate = 0;
4136 if (xRotationSymmetry) {
4137 rcoordinate = sqrt(ypos * ypos + zpos * zpos);
4138 zcoordinate = xpos;
4139 } else if (yRotationSymmetry) {
4140 rcoordinate = sqrt(xpos * xpos + zpos * zpos);
4141 zcoordinate = ypos;
4142 } else if (zRotationSymmetry) {
4143 rcoordinate = sqrt(xpos * xpos + ypos * ypos);
4144 zcoordinate = zpos;
4145 }
4146
4147 if (xRotationSymmetry || yRotationSymmetry || zRotationSymmetry) {
4148 xpos = rcoordinate;
4149 ypos = zcoordinate;
4150 zpos = 0;
4151 }
4152}
cos
DoubleAc cos(const DoubleAc &f)
Definition: DoubleAc.cpp:431
sin
DoubleAc sin(const DoubleAc &f)
Definition: DoubleAc.cpp:383

Referenced by Garfield::ComponentCST::Coordinate2Index(), Garfield::ComponentAnsys121::ElectricField(), Garfield::ComponentAnsys123::ElectricField(), Garfield::ComponentElmer::ElectricField(), Garfield::ComponentAnsys121::GetMedium(), Garfield::ComponentAnsys123::GetMedium(), Garfield::ComponentElmer::GetMedium(), Garfield::ComponentAnsys121::WeightingField(), Garfield::ComponentAnsys123::WeightingField(), Garfield::ComponentElmer::WeightingField(), Garfield::ComponentAnsys121::WeightingPotential(), Garfield::ComponentAnsys123::WeightingPotential(), and Garfield::ComponentElmer::WeightingPotential().

◆ NotDriftMedium()

void Garfield::ComponentFieldMap::NotDriftMedium ( int  imat)

Definition at line 92 of file ComponentFieldMap.cc.

92 {
93
94 // Do not proceed if not properly initialised.
95 if (!ready) {
96 std::cerr << m_className << "::NotDriftMedium:\n";
97 std::cerr << " Field map not yet initialised.\n";
98 std::cerr << " Drift medium cannot be selected.\n";
99 return;
100 }
101
102 // Check value
103 if (imat < 0 || imat >= nMaterials) {
104 std::cerr << m_className << "::NotDriftMedium:\n";
105 std::cerr << " Material index " << imat << " is out of range.\n";
106 return;
107 }
108
109 // Make drift medium
110 materials[imat].driftmedium = false;
111}

◆ PrintMaterials()

void Garfield::ComponentFieldMap::PrintMaterials ( )

Definition at line 37 of file ComponentFieldMap.cc.

37 {
38
39 // Do not proceed if not properly initialised.
40 if (!ready) {
41 std::cerr << m_className << "::PrintMaterials:\n";
42 std::cerr << " Field map not yet initialised.\n";
43 return;
44 }
45
46 if (nMaterials < 0) {
47 std::cerr << m_className << "::PrintMaterials:\n";
48 std::cerr << " No materials are currently defined.\n";
49 return;
50 }
51
52 std::cout << m_className << "::PrintMaterials:\n";
53 std::cout << " Currently " << nMaterials << " materials are defined.\n";
54 std::cout << " Index Permittivity Resistivity Notes\n";
55 for (int i = 0; i < nMaterials; ++i) {
56 printf(" %5d %12g %12g", i, materials[i].eps, materials[i].ohm);
57 if (materials[i].medium != 0) {
58 std::string name = materials[i].medium->GetName();
59 std::cout << " " << name;
60 if (materials[i].medium->IsDriftable()) std::cout << ", drift medium";
61 if (materials[i].medium->IsIonisable()) std::cout << ", ionisable";
62 }
63 if (materials[i].driftmedium) {
64 std::cout << " (drift medium)\n";
65 } else {
66 std::cout << "\n";
67 }
68 }
69}
Garfield::Medium::IsIonisable
bool IsIonisable() const
Definition: Medium.hh:59
Garfield::Medium::IsDriftable
bool IsDriftable() const
Definition: Medium.hh:57

Referenced by Garfield::ComponentCST::Initialise(), Garfield::ComponentAnsys121::Initialise(), and Garfield::ComponentAnsys123::Initialise().

◆ PrintRange()

void Garfield::ComponentFieldMap::PrintRange ( )

Definition at line 3950 of file ComponentFieldMap.cc.

3950 {
3951
3952 std::cout << m_className << "::PrintRange:\n";
3953 std::cout << " Dimensions of the elementary block\n";
3954 printf(" %15g < x < %-15g cm,\n", mapxmin, mapxmax);
3955 printf(" %15g < y < %-15g cm,\n", mapymin, mapymax);
3956 printf(" %15g < z < %-15g cm,\n", mapzmin, mapzmax);
3957 printf(" %15g < V < %-15g V.\n", mapvmin, mapvmax);
3958
3959 std::cout << " Periodicities\n";
3960
3961 std::cout << " x:";
3962 if (xPeriodic) {
3963 std::cout << " simple with length " << cellsx << " cm";
3964 }
3965 if (xMirrorPeriodic) {
3966 std::cout << " mirror with length " << cellsx << " cm";
3967 }
3968 if (xAxiallyPeriodic) {
3969 std::cout << " axial " << int(0.5 + mapnxa) << "-fold repetition";
3970 }
3971 if (xRotationSymmetry) std::cout << " rotational symmetry";
3972 if (!(xPeriodic || xMirrorPeriodic || xAxiallyPeriodic || xRotationSymmetry))
3973 std::cout << " none";
3974 std::cout << "\n";
3975
3976 std::cout << " y:";
3977 if (yPeriodic) {
3978 std::cout << " simple with length " << cellsy << " cm";
3979 }
3980 if (yMirrorPeriodic) {
3981 std::cout << " mirror with length " << cellsy << " cm";
3982 }
3983 if (yAxiallyPeriodic) {
3984 std::cout << " axial " << int(0.5 + mapnya) << "-fold repetition";
3985 }
3986 if (yRotationSymmetry) {
3987 std::cout << " rotational symmetry";
3988 }
3989 if (!(yPeriodic || yMirrorPeriodic || yAxiallyPeriodic || yRotationSymmetry))
3990 std::cout << " none";
3991 std::cout << "\n";
3992
3993 std::cout << " z:";
3994 if (zPeriodic) {
3995 std::cout << " simple with length " << cellsz << " cm";
3996 }
3997 if (zMirrorPeriodic) {
3998 std::cout << " mirror with length " << cellsz << " cm";
3999 }
4000 if (zAxiallyPeriodic) {
4001 std::cout << " axial " << int(0.5 + mapnza) << "-fold repetition";
4002 }
4003 if (zRotationSymmetry) {
4004 std::cout << " rotational symmetry";
4005 }
4006 if (!(zPeriodic || zMirrorPeriodic || zAxiallyPeriodic || zRotationSymmetry))
4007 std::cout << " none";
4008 std::cout << "\n";
4009}

Referenced by Garfield::ComponentCST::Coordinate2Index(), SetRange(), and UpdatePeriodicityCommon().

◆ ReadDouble()

double Garfield::ComponentFieldMap::ReadDouble ( char *  token,
double  def,
bool &  error 
)
protected

Definition at line 4240 of file ComponentFieldMap.cc.

4240 {
4241
4242 if (!token) {
4243 error = true;
4244 return def;
4245 }
4246 return atof(token);
4247}

Referenced by Garfield::ComponentCST::Initialise(), Garfield::ComponentAnsys121::Initialise(), Garfield::ComponentAnsys123::Initialise(), Garfield::ComponentElmer::Initialise(), Garfield::ComponentAnsys121::SetWeightingField(), Garfield::ComponentAnsys123::SetWeightingField(), Garfield::ComponentElmer::SetWeightingField(), and Garfield::ComponentCST::SetWeightingField().

◆ ReadInteger()

int Garfield::ComponentFieldMap::ReadInteger ( char *  token,
int  def,
bool &  error 
)
protected

Definition at line 4230 of file ComponentFieldMap.cc.

4230 {
4231
4232 if (!token) {
4233 error = true;
4234 return def;
4235 }
4236
4237 return atoi(token);
4238}

Referenced by Garfield::ComponentCST::Initialise(), Garfield::ComponentAnsys121::Initialise(), Garfield::ComponentAnsys123::Initialise(), Garfield::ComponentElmer::Initialise(), Garfield::ComponentAnsys121::SetWeightingField(), Garfield::ComponentAnsys123::SetWeightingField(), and Garfield::ComponentCST::SetWeightingField().

◆ Reset()

void Garfield::ComponentFieldMap::Reset ( )
inlineprotectedvirtual

Implements Garfield::ComponentBase.

Definition at line 158 of file ComponentFieldMap.hh.

158{};

◆ SetMedium()

void Garfield::ComponentFieldMap::SetMedium ( const int  imat,
Medium *  medium 
)

Definition at line 135 of file ComponentFieldMap.cc.

135 {
136
137 if (imat < 0 || imat >= nMaterials) {
138 std::cerr << m_className << "::SetMedium:\n";
139 std::cerr << " Material index " << imat << " is out of range.\n";
140 return;
141 }
142
143 if (!m) {
144 std::cerr << m_className << "::SetMedium:\n";
145 std::cerr << " Medium pointer is null.\n";
146 return;
147 }
148
149 if (debug) {
150 std::string name = m->GetName();
151 std::cout << m_className << "::SetMedium:\n";
152 std::cout << " Associated material " << imat << " with medium " << name
153 << ".\n";
154 }
155
156 materials[imat].medium = m;
157}

◆ SetRange()

void Garfield::ComponentFieldMap::SetRange ( )
virtual

Reimplemented in Garfield::ComponentCST.

Definition at line 3835 of file ComponentFieldMap.cc.

3835 {
3836
3837 // Initial values
3838 mapxmin = mapymin = mapzmin = 0.;
3839 mapxmax = mapymax = mapzmax = 0.;
3840 mapxamin = mapyamin = mapzamin = 0.;
3841 mapxamax = mapyamax = mapzamax = 0.;
3842 mapvmin = mapvmax = 0.;
3843 setangx = setangy = setangz = false;
3844
3845 // Make sure the required data is available.
3846 if (!ready || nNodes < 1) {
3847 std::cerr << m_className << "::SetRange:\n";
3848 std::cerr << " Field map not yet set.\n";
3849 return;
3850 }
3851 if (nNodes < 1) {
3852 std::cerr << m_className << "::SetRange:\n";
3853 std::cerr << " Number of nodes < 1.\n";
3854 return;
3855 }
3856
3857 // Loop over the nodes.
3858 mapxmin = mapxmax = nodes[0].x;
3859 mapymin = mapymax = nodes[0].y;
3860 mapzmin = mapzmax = nodes[0].z;
3861 mapvmin = mapvmax = nodes[0].v;
3862
3863 double ang;
3864 for (int i = 1; i < nNodes; i++) {
3865 if (mapxmin > nodes[i].x) mapxmin = nodes[i].x;
3866 if (mapxmax < nodes[i].x) mapxmax = nodes[i].x;
3867 if (mapymin > nodes[i].y) mapymin = nodes[i].y;
3868 if (mapymax < nodes[i].y) mapymax = nodes[i].y;
3869 if (mapzmin > nodes[i].z) mapzmin = nodes[i].z;
3870 if (mapzmax < nodes[i].z) mapzmax = nodes[i].z;
3871 if (mapvmin > nodes[i].v) mapvmin = nodes[i].v;
3872 if (mapvmax < nodes[i].v) mapvmax = nodes[i].v;
3873
3874 if (nodes[i].y != 0 || nodes[i].z != 0) {
3875 ang = atan2(nodes[i].z, nodes[i].y);
3876 if (setangx) {
3877 if (ang < mapxamin) mapxamin = ang;
3878 if (ang > mapxamax) mapxamax = ang;
3879 } else {
3880 mapxamin = mapxamax = ang;
3881 setangx = true;
3882 }
3883 }
3884
3885 if (nodes[i].z != 0 || nodes[i].x != 0) {
3886 ang = atan2(nodes[i].x, nodes[i].z);
3887 if (setangy) {
3888 if (ang < mapyamin) mapyamin = ang;
3889 if (ang > mapyamax) mapyamax = ang;
3890 } else {
3891 mapyamin = mapyamax = ang;
3892 setangy = true;
3893 }
3894 }
3895
3896 if (nodes[i].x != 0 || nodes[i].y != 0) {
3897 ang = atan2(nodes[i].y, nodes[i].x);
3898 if (setangz) {
3899 if (ang < mapzamin) mapzamin = ang;
3900 if (ang > mapzamax) mapzamax = ang;
3901 } else {
3902 mapzamin = mapzamax = ang;
3903 setangz = true;
3904 }
3905 }
3906 }
3907
3908 // Fix the angular ranges.
3909 if (mapxamax - mapxamin > Pi) {
3910 double aux = mapxamin;
3911 mapxamin = mapxamax;
3912 mapxamax = aux + TwoPi;
3913 }
3914
3915 if (mapyamax - mapyamin > Pi) {
3916 double aux = mapyamin;
3917 mapyamin = mapyamax;
3918 mapyamax = aux + TwoPi;
3919 }
3920
3921 if (mapzamax - mapzamin > Pi) {
3922 double aux = mapzamin;
3923 mapzamin = mapzamax;
3924 mapzamax = aux + TwoPi;
3925 }
3926
3927 // Set the periodicity length (maybe not needed).
3928 mapsx = fabs(mapxmax - mapxmin);
3929 mapsy = fabs(mapymax - mapymin);
3930 mapsz = fabs(mapzmax - mapzmin);
3931
3932 // Set provisional cell dimensions.
3933 xMinBoundingBox = mapxmin;
3934 xMaxBoundingBox = mapxmax;
3935 yMinBoundingBox = mapymin;
3936 yMaxBoundingBox = mapymax;
3937 if (is3d) {
3938 zMinBoundingBox = mapzmin;
3939 zMaxBoundingBox = mapzmax;
3940 } else {
3941 mapzmin = zMinBoundingBox;
3942 mapzmax = zMaxBoundingBox;
3943 }
3944 hasBoundingBox = true;
3945
3946 // Display the range if requested.
3947 if (debug) PrintRange();
3948}

Referenced by Garfield::ComponentAnsys121::Initialise(), Garfield::ComponentAnsys123::Initialise(), and Garfield::ComponentElmer::Initialise().

◆ UnmapFields()

void Garfield::ComponentFieldMap::UnmapFields ( double &  ex,
double &  ey,
double &  ez,
double &  xpos,
double &  ypos,
double &  zpos,
bool &  xmirrored,
bool &  ymirrored,
bool &  zmirrored,
double &  rcoordinate,
double &  rotation 
)
protected

Definition at line 4154 of file ComponentFieldMap.cc.

4158 {
4159
4160 // Apply mirror imaging.
4161 if (xmirrored) ex = -ex;
4162 if (ymirrored) ey = -ey;
4163 if (zmirrored) ez = -ez;
4164
4165 // Rotate the field.
4166 double er, theta;
4167 if (xAxiallyPeriodic) {
4168 er = sqrt(ey * ey + ez * ez);
4169 theta = atan2(ez, ey);
4170 theta += rotation;
4171 ey = er * cos(theta);
4172 ez = er * sin(theta);
4173 }
4174 if (yAxiallyPeriodic) {
4175 er = sqrt(ez * ez + ex * ex);
4176 theta = atan2(ex, ez);
4177 theta += rotation;
4178 ez = er * cos(theta);
4179 ex = er * sin(theta);
4180 }
4181 if (zAxiallyPeriodic) {
4182 er = sqrt(ex * ex + ey * ey);
4183 theta = atan2(ey, ex);
4184 theta += rotation;
4185 ex = er * cos(theta);
4186 ey = er * sin(theta);
4187 }
4188
4189 // Take care of symmetry.
4190 double eaxis;
4191 er = ex;
4192 eaxis = ey;
4193
4194 // Rotational symmetry
4195 if (xRotationSymmetry) {
4196 if (rcoordinate <= 0) {
4197 ex = eaxis;
4198 ey = 0;
4199 ez = 0;
4200 } else {
4201 ex = eaxis;
4202 ey = er * ypos / rcoordinate;
4203 ez = er * zpos / rcoordinate;
4204 }
4205 }
4206 if (yRotationSymmetry) {
4207 if (rcoordinate <= 0) {
4208 ex = 0;
4209 ey = eaxis;
4210 ez = 0;
4211 } else {
4212 ex = er * xpos / rcoordinate;
4213 ey = eaxis;
4214 ez = er * zpos / rcoordinate;
4215 }
4216 }
4217 if (zRotationSymmetry) {
4218 if (rcoordinate <= 0) {
4219 ex = 0;
4220 ey = 0;
4221 ez = eaxis;
4222 } else {
4223 ex = er * xpos / rcoordinate;
4224 ey = er * ypos / rcoordinate;
4225 ez = eaxis;
4226 }
4227 }
4228}

Referenced by Garfield::ComponentAnsys121::ElectricField(), Garfield::ComponentAnsys123::ElectricField(), Garfield::ComponentElmer::ElectricField(), Garfield::ComponentAnsys121::WeightingField(), Garfield::ComponentAnsys123::WeightingField(), and Garfield::ComponentElmer::WeightingField().

◆ UpdatePeriodicity()

virtual void Garfield::ComponentFieldMap::UpdatePeriodicity ( )
protectedpure virtual

Implements Garfield::ComponentBase.

Implemented in Garfield::ComponentAnsys121, Garfield::ComponentAnsys123, Garfield::ComponentCST, and Garfield::ComponentElmer.

◆ UpdatePeriodicity2d()

void Garfield::ComponentFieldMap::UpdatePeriodicity2d ( )
protected

Definition at line 3805 of file ComponentFieldMap.cc.

3805 {
3806
3807 // Check the required data is available.
3808 if (!ready) {
3809 std::cerr << m_className << "::UpdatePeriodicity2d:\n";
3810 std::cerr << " No valid field map available.\n";
3811 return;
3812 }
3813
3814 // No z-periodicity in 2d
3815 if (zPeriodic || zMirrorPeriodic) {
3816 std::cerr << m_className << "::UpdatePeriodicity2d:\n";
3817 std::cerr << " Simple or mirror periodicity along z\n";
3818 std::cerr << " requested for a 2d map; reset.\n";
3819 zPeriodic = false;
3820 zMirrorPeriodic = false;
3821 warning = true;
3822 }
3823
3824 // Only z-axial periodicity in 2d maps
3825 if (xAxiallyPeriodic || yAxiallyPeriodic) {
3826 std::cerr << m_className << "::UpdatePeriodicity2d:\n";
3827 std::cerr << " Axial symmetry has been requested \n";
3828 std::cerr << " around x or y for a 2D map; reset.\n";
3829 xAxiallyPeriodic = false;
3830 yAxiallyPeriodic = false;
3831 warning = true;
3832 }
3833}

Referenced by Garfield::ComponentAnsys121::UpdatePeriodicity(), and Garfield::ComponentCST::UpdatePeriodicity().

◆ UpdatePeriodicityCommon()

void Garfield::ComponentFieldMap::UpdatePeriodicityCommon ( )
protected

Definition at line 3604 of file ComponentFieldMap.cc.

3604 {
3605
3606 // Check the required data is available.
3607 if (!ready) {
3608 std::cerr << m_className << "::UpdatePeriodicityCommon:\n";
3609 std::cerr << " No valid field map available.\n";
3610 return;
3611 }
3612
3613 // No regular and mirror periodicity at the same time.
3614 if (xPeriodic && xMirrorPeriodic) {
3615 std::cerr << m_className << "::UpdatePeriodicityCommon:\n";
3616 std::cerr << " Both simple and mirror periodicity\n";
3617 std::cerr << " along x requested; reset.\n";
3618 xPeriodic = false;
3619 xMirrorPeriodic = false;
3620 warning = true;
3621 }
3622 if (yPeriodic && yMirrorPeriodic) {
3623 std::cerr << m_className << "::UpdatePeriodicityCommon:\n";
3624 std::cerr << " Both simple and mirror periodicity\n";
3625 std::cerr << " along y requested; reset.\n";
3626 yPeriodic = false;
3627 yMirrorPeriodic = false;
3628 warning = true;
3629 }
3630 if (zPeriodic && zMirrorPeriodic) {
3631 std::cerr << m_className << "::UpdatePeriodicityCommon:\n";
3632 std::cerr << " Both simple and mirror periodicity\n";
3633 std::cerr << " along z requested; reset.\n";
3634 zPeriodic = false;
3635 zMirrorPeriodic = false;
3636 warning = true;
3637 }
3638
3639 // In case of axial periodicity,
3640 // the range must be an integral part of 2 pi.
3641 if (xAxiallyPeriodic) {
3642 if (mapxamin >= mapxamax) {
3643 mapnxa = 0;
3644 } else {
3645 mapnxa = TwoPi / (mapxamax - mapxamin);
3646 }
3647 if (fabs(mapnxa - int(0.5 + mapnxa)) > 0.001 || mapnxa < 1.5) {
3648 std::cerr << m_className << "::UpdatePeriodicityCommon:\n";
3649 std::cerr << " X-axial symmetry has been requested but the map\n";
3650 std::cerr << " does not cover an integral fraction of 2 pi; reset.\n";
3651 xAxiallyPeriodic = false;
3652 warning = true;
3653 }
3654 }
3655
3656 if (yAxiallyPeriodic) {
3657 if (mapyamin >= mapyamax) {
3658 mapnya = 0;
3659 } else {
3660 mapnya = TwoPi / (mapyamax - mapyamin);
3661 }
3662 if (fabs(mapnya - int(0.5 + mapnya)) > 0.001 || mapnya < 1.5) {
3663 std::cerr << m_className << "::UpdatePeriodicityCommon:\n";
3664 std::cerr << " Y-axial symmetry has been requested but the map\n";
3665 std::cerr << " does not cover an integral fraction of 2 pi; reset.\n";
3666 yAxiallyPeriodic = false;
3667 warning = true;
3668 }
3669 }
3670
3671 if (zAxiallyPeriodic) {
3672 if (mapzamin >= mapzamax) {
3673 mapnza = 0;
3674 } else {
3675 mapnza = TwoPi / (mapzamax - mapzamin);
3676 }
3677 if (fabs(mapnza - int(0.5 + mapnza)) > 0.001 || mapnza < 1.5) {
3678 std::cerr << m_className << "::UpdatePeriodicityCommon:\n";
3679 std::cerr << " Z-axial symmetry has been requested but the map\n";
3680 std::cerr << " does not cover an integral fraction of 2 pi; reset.\n";
3681 zAxiallyPeriodic = false;
3682 warning = true;
3683 }
3684 }
3685
3686 // Not more than 1 rotational symmetry
3687 if ((xRotationSymmetry && yRotationSymmetry) ||
3688 (xRotationSymmetry && zRotationSymmetry) ||
3689 (yRotationSymmetry && zRotationSymmetry)) {
3690 std::cerr << m_className << "::UpdatePeriodicityCommon:\n";
3691 std::cerr << " Only 1 rotational symmetry allowed; reset.\n";
3692 xRotationSymmetry = false;
3693 yRotationSymmetry = false;
3694 zRotationSymmetry = false;
3695 warning = true;
3696 }
3697
3698 // No rotational symmetry as well as axial periodicity
3699 if ((xRotationSymmetry || yRotationSymmetry || zRotationSymmetry) &&
3700 (xAxiallyPeriodic || yAxiallyPeriodic || zAxiallyPeriodic)) {
3701 std::cerr << m_className << "::UpdatePeriodicityCommon:\n";
3702 std::cerr << " Not allowed to combine rotational symmetry\n";
3703 std::cerr << " and axial periodicity; reset.\n";
3704 xAxiallyPeriodic = false;
3705 yAxiallyPeriodic = false;
3706 zAxiallyPeriodic = false;
3707 xRotationSymmetry = false;
3708 yRotationSymmetry = false;
3709 zRotationSymmetry = false;
3710 warning = true;
3711 }
3712
3713 // In case of rotational symmetry, the x-range should not straddle 0.
3714 if (xRotationSymmetry || yRotationSymmetry || zRotationSymmetry) {
3715 if (mapxmin * mapxmax < 0) {
3716 std::cerr << m_className << "::UpdatePeriodicityCommon:\n";
3717 std::cerr << " Rotational symmetry requested, \n";
3718 std::cerr << " but x-range straddles 0; reset.\n";
3719 xRotationSymmetry = false;
3720 yRotationSymmetry = false;
3721 zRotationSymmetry = false;
3722 warning = true;
3723 }
3724 }
3725
3726 // Recompute the cell ranges.
3727 xMinBoundingBox = mapxmin;
3728 xMaxBoundingBox = mapxmax;
3729 yMinBoundingBox = mapymin;
3730 yMaxBoundingBox = mapymax;
3731 zMinBoundingBox = mapzmin;
3732 zMaxBoundingBox = mapzmax;
3733 cellsx = fabs(mapxmax - mapxmin);
3734 cellsy = fabs(mapymax - mapymin);
3735 cellsz = fabs(mapzmax - mapzmin);
3736 if (xRotationSymmetry) {
3737 xMinBoundingBox = mapymin;
3738 xMaxBoundingBox = mapymax;
3739 yMinBoundingBox = -std::max(fabs(mapxmin), fabs(mapxmax));
3740 yMaxBoundingBox = +std::max(fabs(mapxmin), fabs(mapxmax));
3741 zMinBoundingBox = -std::max(fabs(mapxmin), fabs(mapxmax));
3742 zMaxBoundingBox = +std::max(fabs(mapxmin), fabs(mapxmax));
3743 } else if (yRotationSymmetry) {
3744 xMinBoundingBox = -std::max(fabs(mapxmin), fabs(mapxmax));
3745 xMaxBoundingBox = +std::max(fabs(mapxmin), fabs(mapxmax));
3746 yMinBoundingBox = mapymin;
3747 yMaxBoundingBox = mapymax;
3748 zMinBoundingBox = -std::max(fabs(mapxmin), fabs(mapxmax));
3749 zMaxBoundingBox = +std::max(fabs(mapxmin), fabs(mapxmax));
3750 } else if (zRotationSymmetry) {
3751 xMinBoundingBox = -std::max(fabs(mapxmin), fabs(mapxmax));
3752 xMaxBoundingBox = +std::max(fabs(mapxmin), fabs(mapxmax));
3753 yMinBoundingBox = -std::max(fabs(mapxmin), fabs(mapxmax));
3754 yMaxBoundingBox = +std::max(fabs(mapxmin), fabs(mapxmax));
3755 zMinBoundingBox = mapymin;
3756 zMaxBoundingBox = mapymax;
3757 }
3758
3759 if (xAxiallyPeriodic) {
3760 yMinBoundingBox = -std::max(std::max(fabs(mapymin), fabs(mapymax)),
3761 std::max(fabs(mapzmin), fabs(mapzmax)));
3762 yMaxBoundingBox = +std::max(std::max(fabs(mapymin), fabs(mapymax)),
3763 std::max(fabs(mapzmin), fabs(mapzmax)));
3764 zMinBoundingBox = -std::max(std::max(fabs(mapymin), fabs(mapymax)),
3765 std::max(fabs(mapzmin), fabs(mapzmax)));
3766 zMaxBoundingBox = +std::max(std::max(fabs(mapymin), fabs(mapymax)),
3767 std::max(fabs(mapzmin), fabs(mapzmax)));
3768 } else if (yAxiallyPeriodic) {
3769 xMinBoundingBox = -std::max(std::max(fabs(mapxmin), fabs(mapxmax)),
3770 std::max(fabs(mapzmin), fabs(mapzmax)));
3771 xMaxBoundingBox = +std::max(std::max(fabs(mapxmin), fabs(mapxmax)),
3772 std::max(fabs(mapzmin), fabs(mapzmax)));
3773 zMinBoundingBox = -std::max(std::max(fabs(mapxmin), fabs(mapxmax)),
3774 std::max(fabs(mapzmin), fabs(mapzmax)));
3775 zMaxBoundingBox = +std::max(std::max(fabs(mapxmin), fabs(mapxmax)),
3776 std::max(fabs(mapzmin), fabs(mapzmax)));
3777 } else if (zAxiallyPeriodic) {
3778 xMinBoundingBox = -std::max(std::max(fabs(mapxmin), fabs(mapxmax)),
3779 std::max(fabs(mapymin), fabs(mapymax)));
3780 xMaxBoundingBox = +std::max(std::max(fabs(mapxmin), fabs(mapxmax)),
3781 std::max(fabs(mapymin), fabs(mapymax)));
3782 yMinBoundingBox = -std::max(std::max(fabs(mapxmin), fabs(mapxmax)),
3783 std::max(fabs(mapymin), fabs(mapymax)));
3784 yMaxBoundingBox = +std::max(std::max(fabs(mapxmin), fabs(mapxmax)),
3785 std::max(fabs(mapymin), fabs(mapymax)));
3786 }
3787
3788 if (xPeriodic || xMirrorPeriodic) {
3789 xMinBoundingBox = -INFINITY;
3790 xMaxBoundingBox = +INFINITY;
3791 }
3792 if (yPeriodic || yMirrorPeriodic) {
3793 yMinBoundingBox = -INFINITY;
3794 yMaxBoundingBox = +INFINITY;
3795 }
3796 if (zPeriodic || zMirrorPeriodic) {
3797 zMinBoundingBox = -INFINITY;
3798 zMaxBoundingBox = +INFINITY;
3799 }
3800
3801 // Display the range if requested.
3802 if (debug) PrintRange();
3803}

Referenced by Garfield::ComponentAnsys121::UpdatePeriodicity(), Garfield::ComponentAnsys123::UpdatePeriodicity(), Garfield::ComponentCST::UpdatePeriodicity(), and Garfield::ComponentElmer::UpdatePeriodicity().

◆ WeightingField()

virtual void Garfield::ComponentFieldMap::WeightingField ( const double  x,
const double  y,
const double  z,
double &  wx,
double &  wy,
double &  wz,
const std::string  label 
)
pure virtual

Reimplemented from Garfield::ComponentBase.

Implemented in Garfield::ComponentAnsys121, Garfield::ComponentAnsys123, Garfield::ComponentCST, and Garfield::ComponentElmer.

◆ WeightingPotential()

virtual double Garfield::ComponentFieldMap::WeightingPotential ( const double  x,
const double  y,
const double  z,
const std::string  label 
)
pure virtual

Reimplemented from Garfield::ComponentBase.

Implemented in Garfield::ComponentAnsys121, Garfield::ComponentAnsys123, Garfield::ComponentCST, and Garfield::ComponentElmer.

Friends And Related Function Documentation

◆ ViewFEMesh

friend class ViewFEMesh
friend

Definition at line 80 of file ComponentFieldMap.hh.

Member Data Documentation

◆ cacheElemBoundingBoxes

bool Garfield::ComponentFieldMap::cacheElemBoundingBoxes
protected

Definition at line 99 of file ComponentFieldMap.hh.

Referenced by FindElement13(), and FindElement5().

◆ cellsx

double Garfield::ComponentFieldMap::cellsx
protected

Definition at line 145 of file ComponentFieldMap.hh.

Referenced by PrintRange(), and UpdatePeriodicityCommon().

◆ cellsy

double Garfield::ComponentFieldMap::cellsy
protected

Definition at line 145 of file ComponentFieldMap.hh.

Referenced by PrintRange(), and UpdatePeriodicityCommon().

◆ cellsz

double Garfield::ComponentFieldMap::cellsz
protected

Definition at line 145 of file ComponentFieldMap.hh.

Referenced by PrintRange(), and UpdatePeriodicityCommon().

◆ checkMultipleElement

bool Garfield::ComponentFieldMap::checkMultipleElement
protected

Definition at line 152 of file ComponentFieldMap.hh.

Referenced by DisableCheckMapIndices(), EnableCheckMapIndices(), FindElement13(), and FindElement5().

◆ deleteBackground

bool Garfield::ComponentFieldMap::deleteBackground
protected

Definition at line 150 of file ComponentFieldMap.hh.

Referenced by Garfield::ComponentCST::ComponentCST(), DisableDeleteBackgroundElements(), EnableDeleteBackgroundElements(), Garfield::ComponentAnsys121::Initialise(), and Garfield::ComponentAnsys123::Initialise().

◆ elements

std::vector<element> Garfield::ComponentFieldMap::elements
protected

Definition at line 96 of file ComponentFieldMap.hh.

Referenced by CalculateElementBoundingBoxes(), ComponentFieldMap(), Coordinates12(), Coordinates13(), Coordinates3(), Coordinates4(), Coordinates5(), CoordinatesCube(), Garfield::ComponentAnsys121::ElectricField(), Garfield::ComponentAnsys123::ElectricField(), Garfield::ComponentElmer::ElectricField(), FindElement13(), FindElement5(), FindElementCube(), Garfield::ComponentAnsys121::GetAspectRatio(), Garfield::ComponentAnsys123::GetAspectRatio(), Garfield::ComponentElmer::GetAspectRatio(), Garfield::ComponentAnsys121::GetElementVolume(), Garfield::ComponentAnsys123::GetElementVolume(), Garfield::ComponentElmer::GetElementVolume(), Garfield::ComponentAnsys121::GetMedium(), Garfield::ComponentAnsys123::GetMedium(), Garfield::ComponentElmer::GetMedium(), Garfield::ComponentAnsys121::Initialise(), Garfield::ComponentAnsys123::Initialise(), Garfield::ComponentElmer::Initialise(), Jacobian13(), Jacobian3(), Jacobian5(), JacobianCube(), Garfield::ComponentAnsys121::WeightingField(), Garfield::ComponentAnsys123::WeightingField(), Garfield::ComponentElmer::WeightingField(), Garfield::ComponentAnsys121::WeightingPotential(), Garfield::ComponentAnsys123::WeightingPotential(), and Garfield::ComponentElmer::WeightingPotential().

◆ hasBoundingBox

bool Garfield::ComponentFieldMap::hasBoundingBox
protected

Definition at line 131 of file ComponentFieldMap.hh.

Referenced by Garfield::ComponentCST::SetRange(), and SetRange().

◆ is3d

bool Garfield::ComponentFieldMap::is3d
protected

Definition at line 83 of file ComponentFieldMap.hh.

Referenced by Garfield::ComponentAnsys121::ComponentAnsys121(), Garfield::ComponentCST::SetRange(), and SetRange().

◆ lastElement

int Garfield::ComponentFieldMap::lastElement
protected

Definition at line 97 of file ComponentFieldMap.hh.

Referenced by EnableCheckMapIndices(), FindElement13(), FindElement5(), and FindElementCube().

◆ mapnxa

double Garfield::ComponentFieldMap::mapnxa
protected

Definition at line 146 of file ComponentFieldMap.hh.

Referenced by PrintRange(), and UpdatePeriodicityCommon().

◆ mapnya

double Garfield::ComponentFieldMap::mapnya
protected

Definition at line 146 of file ComponentFieldMap.hh.

Referenced by PrintRange(), and UpdatePeriodicityCommon().

◆ mapnza

double Garfield::ComponentFieldMap::mapnza
protected

Definition at line 146 of file ComponentFieldMap.hh.

Referenced by PrintRange(), and UpdatePeriodicityCommon().

◆ mapsx

double Garfield::ComponentFieldMap::mapsx
protected

Definition at line 143 of file ComponentFieldMap.hh.

Referenced by Garfield::ComponentCST::SetRange(), and SetRange().

◆ mapsy

double Garfield::ComponentFieldMap::mapsy
protected

Definition at line 143 of file ComponentFieldMap.hh.

Referenced by Garfield::ComponentCST::SetRange(), and SetRange().

◆ mapsz

double Garfield::ComponentFieldMap::mapsz
protected

Definition at line 143 of file ComponentFieldMap.hh.

Referenced by Garfield::ComponentCST::SetRange(), and SetRange().

◆ mapvmax

double Garfield::ComponentFieldMap::mapvmax
protected

Definition at line 140 of file ComponentFieldMap.hh.

Referenced by GetVoltageRange(), PrintRange(), Garfield::ComponentCST::SetRange(), and SetRange().

◆ mapvmin

double Garfield::ComponentFieldMap::mapvmin
protected

Definition at line 140 of file ComponentFieldMap.hh.

Referenced by GetVoltageRange(), PrintRange(), Garfield::ComponentCST::SetRange(), and SetRange().

◆ mapxamax

double Garfield::ComponentFieldMap::mapxamax
protected

Definition at line 139 of file ComponentFieldMap.hh.

Referenced by MapCoordinates(), SetRange(), and UpdatePeriodicityCommon().

◆ mapxamin

double Garfield::ComponentFieldMap::mapxamin
protected

Definition at line 138 of file ComponentFieldMap.hh.

Referenced by MapCoordinates(), SetRange(), and UpdatePeriodicityCommon().

◆ mapxmax

double Garfield::ComponentFieldMap::mapxmax
protected

Definition at line 137 of file ComponentFieldMap.hh.

Referenced by MapCoordinates(), PrintRange(), Garfield::ComponentCST::SetRange(), SetRange(), and UpdatePeriodicityCommon().

◆ mapxmin

double Garfield::ComponentFieldMap::mapxmin
protected

Definition at line 136 of file ComponentFieldMap.hh.

Referenced by MapCoordinates(), PrintRange(), Garfield::ComponentCST::SetRange(), SetRange(), and UpdatePeriodicityCommon().

◆ mapyamax

double Garfield::ComponentFieldMap::mapyamax
protected

Definition at line 139 of file ComponentFieldMap.hh.

Referenced by MapCoordinates(), SetRange(), and UpdatePeriodicityCommon().

◆ mapyamin

double Garfield::ComponentFieldMap::mapyamin
protected

Definition at line 138 of file ComponentFieldMap.hh.

Referenced by MapCoordinates(), SetRange(), and UpdatePeriodicityCommon().

◆ mapymax

double Garfield::ComponentFieldMap::mapymax
protected

Definition at line 137 of file ComponentFieldMap.hh.

Referenced by MapCoordinates(), PrintRange(), Garfield::ComponentCST::SetRange(), SetRange(), and UpdatePeriodicityCommon().

◆ mapymin

double Garfield::ComponentFieldMap::mapymin
protected

Definition at line 136 of file ComponentFieldMap.hh.

Referenced by MapCoordinates(), PrintRange(), Garfield::ComponentCST::SetRange(), SetRange(), and UpdatePeriodicityCommon().

◆ mapzamax

double Garfield::ComponentFieldMap::mapzamax
protected

Definition at line 139 of file ComponentFieldMap.hh.

Referenced by MapCoordinates(), SetRange(), and UpdatePeriodicityCommon().

◆ mapzamin

double Garfield::ComponentFieldMap::mapzamin
protected

Definition at line 138 of file ComponentFieldMap.hh.

Referenced by MapCoordinates(), SetRange(), and UpdatePeriodicityCommon().

◆ mapzmax

double Garfield::ComponentFieldMap::mapzmax
protected

Definition at line 137 of file ComponentFieldMap.hh.

Referenced by MapCoordinates(), PrintRange(), Garfield::ComponentCST::SetRange(), SetRange(), Garfield::ComponentAnsys121::SetRangeZ(), and UpdatePeriodicityCommon().

◆ mapzmin

double Garfield::ComponentFieldMap::mapzmin
protected

Definition at line 136 of file ComponentFieldMap.hh.

Referenced by MapCoordinates(), PrintRange(), Garfield::ComponentCST::SetRange(), SetRange(), Garfield::ComponentAnsys121::SetRangeZ(), and UpdatePeriodicityCommon().

◆ materials

std::vector<material> Garfield::ComponentFieldMap::materials
protected

Definition at line 124 of file ComponentFieldMap.hh.

Referenced by ComponentFieldMap(), DriftMedium(), Garfield::ComponentAnsys121::ElectricField(), Garfield::ComponentAnsys123::ElectricField(), Garfield::ComponentElmer::ElectricField(), GetConductivity(), Garfield::ComponentAnsys121::GetMedium(), Garfield::ComponentAnsys123::GetMedium(), Garfield::ComponentCST::GetMedium(), Garfield::ComponentElmer::GetMedium(), GetMedium(), GetPermittivity(), Garfield::ComponentCST::Initialise(), Garfield::ComponentAnsys121::Initialise(), Garfield::ComponentAnsys123::Initialise(), Garfield::ComponentElmer::Initialise(), NotDriftMedium(), PrintMaterials(), SetMedium(), and Garfield::ComponentCST::WeightingField().

◆ nElements

int Garfield::ComponentFieldMap::nElements
protected

Definition at line 86 of file ComponentFieldMap.hh.

Referenced by CalculateElementBoundingBoxes(), FindElement13(), FindElement5(), FindElementCube(), Garfield::ComponentAnsys121::GetAspectRatio(), Garfield::ComponentAnsys123::GetAspectRatio(), Garfield::ComponentCST::GetAspectRatio(), Garfield::ComponentElmer::GetAspectRatio(), GetElement(), Garfield::ComponentAnsys121::GetElementVolume(), Garfield::ComponentAnsys123::GetElementVolume(), Garfield::ComponentCST::GetElementVolume(), Garfield::ComponentElmer::GetElementVolume(), GetNumberOfElements(), Garfield::ComponentCST::Initialise(), Garfield::ComponentAnsys121::Initialise(), Garfield::ComponentAnsys123::Initialise(), Garfield::ComponentElmer::Initialise(), Jacobian13(), Jacobian3(), Jacobian5(), and JacobianCube().

◆ nMaterials

int Garfield::ComponentFieldMap::nMaterials
protected

Definition at line 114 of file ComponentFieldMap.hh.

Referenced by DriftMedium(), GetConductivity(), Garfield::ComponentAnsys121::GetMedium(), Garfield::ComponentAnsys123::GetMedium(), Garfield::ComponentElmer::GetMedium(), GetMedium(), GetNumberOfMaterials(), GetNumberOfMedia(), GetPermittivity(), Garfield::ComponentCST::Initialise(), Garfield::ComponentAnsys121::Initialise(), Garfield::ComponentAnsys123::Initialise(), Garfield::ComponentElmer::Initialise(), NotDriftMedium(), PrintMaterials(), SetMedium(), and Garfield::ComponentCST::SetWeightingField().

◆ nNodes

int Garfield::ComponentFieldMap::nNodes
protected

Definition at line 102 of file ComponentFieldMap.hh.

Referenced by Garfield::ComponentCST::Initialise(), Garfield::ComponentAnsys121::Initialise(), Garfield::ComponentAnsys123::Initialise(), Garfield::ComponentElmer::Initialise(), SetRange(), Garfield::ComponentAnsys121::SetWeightingField(), Garfield::ComponentAnsys123::SetWeightingField(), Garfield::ComponentElmer::SetWeightingField(), and Garfield::ComponentCST::SetWeightingField().

◆ nodes

std::vector<node> Garfield::ComponentFieldMap::nodes
protected

Definition at line 111 of file ComponentFieldMap.hh.

Referenced by CalculateElementBoundingBoxes(), ComponentFieldMap(), Coordinates12(), Coordinates13(), Coordinates3(), Coordinates4(), Coordinates5(), CoordinatesCube(), Garfield::ComponentAnsys121::ElectricField(), Garfield::ComponentAnsys123::ElectricField(), Garfield::ComponentElmer::ElectricField(), FindElement13(), FindElement5(), FindElementCube(), Garfield::ComponentAnsys121::GetAspectRatio(), Garfield::ComponentAnsys123::GetAspectRatio(), Garfield::ComponentElmer::GetAspectRatio(), Garfield::ComponentAnsys121::GetElementVolume(), Garfield::ComponentAnsys123::GetElementVolume(), Garfield::ComponentElmer::GetElementVolume(), Garfield::ComponentAnsys121::GetMedium(), Garfield::ComponentAnsys123::GetMedium(), Garfield::ComponentElmer::GetMedium(), Garfield::ComponentAnsys121::Initialise(), Garfield::ComponentAnsys123::Initialise(), Garfield::ComponentElmer::Initialise(), Jacobian13(), Jacobian3(), Jacobian5(), JacobianCube(), SetRange(), Garfield::ComponentAnsys121::SetWeightingField(), Garfield::ComponentAnsys123::SetWeightingField(), Garfield::ComponentElmer::SetWeightingField(), Garfield::ComponentAnsys121::WeightingField(), Garfield::ComponentAnsys123::WeightingField(), Garfield::ComponentElmer::WeightingField(), Garfield::ComponentAnsys121::WeightingPotential(), Garfield::ComponentAnsys123::WeightingPotential(), and Garfield::ComponentElmer::WeightingPotential().

◆ nWeightingFields

int Garfield::ComponentFieldMap::nWeightingFields
protected

Definition at line 126 of file ComponentFieldMap.hh.

Referenced by Garfield::ComponentAnsys121::Initialise(), Garfield::ComponentAnsys123::Initialise(), Garfield::ComponentElmer::Initialise(), Garfield::ComponentAnsys121::SetWeightingField(), Garfield::ComponentAnsys123::SetWeightingField(), Garfield::ComponentElmer::SetWeightingField(), Garfield::ComponentAnsys121::WeightingField(), Garfield::ComponentAnsys123::WeightingField(), Garfield::ComponentElmer::WeightingField(), Garfield::ComponentAnsys121::WeightingPotential(), Garfield::ComponentAnsys123::WeightingPotential(), and Garfield::ComponentElmer::WeightingPotential().

◆ setangx

bool Garfield::ComponentFieldMap::setangx
protected

Definition at line 142 of file ComponentFieldMap.hh.

Referenced by SetRange().

◆ setangy

bool Garfield::ComponentFieldMap::setangy
protected

Definition at line 142 of file ComponentFieldMap.hh.

Referenced by SetRange().

◆ setangz

bool Garfield::ComponentFieldMap::setangz
protected

Definition at line 142 of file ComponentFieldMap.hh.

Referenced by SetRange().

◆ warning

bool Garfield::ComponentFieldMap::warning
protected

Definition at line 155 of file ComponentFieldMap.hh.

Referenced by Garfield::ComponentAnsys121::ElectricField(), Garfield::ComponentAnsys123::ElectricField(), Garfield::ComponentElmer::ElectricField(), Garfield::ComponentAnsys121::GetMedium(), Garfield::ComponentAnsys123::GetMedium(), Garfield::ComponentElmer::GetMedium(), UpdatePeriodicity2d(), UpdatePeriodicityCommon(), Garfield::ComponentAnsys121::WeightingField(), Garfield::ComponentAnsys123::WeightingField(), Garfield::ComponentCST::WeightingField(), Garfield::ComponentElmer::WeightingField(), Garfield::ComponentAnsys121::WeightingPotential(), Garfield::ComponentAnsys123::WeightingPotential(), and Garfield::ComponentElmer::WeightingPotential().

◆ wfields

std::vector<std::string> Garfield::ComponentFieldMap::wfields
protected

Definition at line 127 of file ComponentFieldMap.hh.

Referenced by ComponentFieldMap(), Garfield::ComponentAnsys121::Initialise(), Garfield::ComponentAnsys123::Initialise(), Garfield::ComponentElmer::Initialise(), Garfield::ComponentAnsys121::SetWeightingField(), Garfield::ComponentAnsys123::SetWeightingField(), Garfield::ComponentElmer::SetWeightingField(), Garfield::ComponentCST::SetWeightingField(), Garfield::ComponentAnsys121::WeightingField(), Garfield::ComponentAnsys123::WeightingField(), Garfield::ComponentElmer::WeightingField(), Garfield::ComponentAnsys121::WeightingPotential(), Garfield::ComponentAnsys123::WeightingPotential(), and Garfield::ComponentElmer::WeightingPotential().

◆ wfieldsOk

std::vector<bool> Garfield::ComponentFieldMap::wfieldsOk
protected

Definition at line 128 of file ComponentFieldMap.hh.

Referenced by ComponentFieldMap(), Garfield::ComponentAnsys121::Initialise(), Garfield::ComponentAnsys123::Initialise(), Garfield::ComponentElmer::Initialise(), Garfield::ComponentAnsys121::SetWeightingField(), Garfield::ComponentAnsys123::SetWeightingField(), Garfield::ComponentElmer::SetWeightingField(), Garfield::ComponentCST::SetWeightingField(), Garfield::ComponentAnsys121::WeightingField(), Garfield::ComponentAnsys123::WeightingField(), Garfield::ComponentCST::WeightingField(), Garfield::ComponentElmer::WeightingField(), Garfield::ComponentAnsys121::WeightingPotential(), Garfield::ComponentAnsys123::WeightingPotential(), Garfield::ComponentCST::WeightingPotential(), and Garfield::ComponentElmer::WeightingPotential().

◆ xMaxBoundingBox

double Garfield::ComponentFieldMap::xMaxBoundingBox
protected

Definition at line 133 of file ComponentFieldMap.hh.

Referenced by GetBoundingBox(), Garfield::ComponentAnsys121::IsInBoundingBox(), Garfield::ComponentAnsys123::IsInBoundingBox(), Garfield::ComponentCST::IsInBoundingBox(), Garfield::ComponentElmer::IsInBoundingBox(), Garfield::ComponentCST::SetRange(), SetRange(), and UpdatePeriodicityCommon().

◆ xMinBoundingBox

double Garfield::ComponentFieldMap::xMinBoundingBox
protected

Definition at line 132 of file ComponentFieldMap.hh.

Referenced by GetBoundingBox(), Garfield::ComponentAnsys121::IsInBoundingBox(), Garfield::ComponentAnsys123::IsInBoundingBox(), Garfield::ComponentCST::IsInBoundingBox(), Garfield::ComponentElmer::IsInBoundingBox(), IsInBoundingBox(), Garfield::ComponentCST::SetRange(), SetRange(), and UpdatePeriodicityCommon().

◆ yMaxBoundingBox

double Garfield::ComponentFieldMap::yMaxBoundingBox
protected

Definition at line 133 of file ComponentFieldMap.hh.

Referenced by GetBoundingBox(), Garfield::ComponentAnsys121::IsInBoundingBox(), Garfield::ComponentAnsys123::IsInBoundingBox(), Garfield::ComponentCST::IsInBoundingBox(), Garfield::ComponentElmer::IsInBoundingBox(), Garfield::ComponentCST::SetRange(), SetRange(), and UpdatePeriodicityCommon().

◆ yMinBoundingBox

double Garfield::ComponentFieldMap::yMinBoundingBox
protected

Definition at line 132 of file ComponentFieldMap.hh.

Referenced by GetBoundingBox(), Garfield::ComponentAnsys121::IsInBoundingBox(), Garfield::ComponentAnsys123::IsInBoundingBox(), Garfield::ComponentCST::IsInBoundingBox(), Garfield::ComponentElmer::IsInBoundingBox(), IsInBoundingBox(), Garfield::ComponentCST::SetRange(), SetRange(), and UpdatePeriodicityCommon().

◆ zMaxBoundingBox

double Garfield::ComponentFieldMap::zMaxBoundingBox
protected

Definition at line 133 of file ComponentFieldMap.hh.

Referenced by Garfield::ComponentAnsys121::ComponentAnsys121(), Garfield::ComponentCST::ComponentCST(), Garfield::ComponentAnsys121::ElectricField(), GetBoundingBox(), Garfield::ComponentAnsys121::GetMedium(), Garfield::ComponentAnsys121::IsInBoundingBox(), Garfield::ComponentAnsys123::IsInBoundingBox(), Garfield::ComponentCST::IsInBoundingBox(), Garfield::ComponentElmer::IsInBoundingBox(), IsInBoundingBox(), Garfield::ComponentCST::SetRange(), SetRange(), Garfield::ComponentAnsys121::SetRangeZ(), Garfield::ComponentCST::SetRangeZ(), and UpdatePeriodicityCommon().

◆ zMinBoundingBox

double Garfield::ComponentFieldMap::zMinBoundingBox
protected

Definition at line 132 of file ComponentFieldMap.hh.

Referenced by Garfield::ComponentAnsys121::ComponentAnsys121(), Garfield::ComponentCST::ComponentCST(), GetBoundingBox(), Garfield::ComponentAnsys121::IsInBoundingBox(), Garfield::ComponentAnsys123::IsInBoundingBox(), Garfield::ComponentCST::IsInBoundingBox(), Garfield::ComponentElmer::IsInBoundingBox(), IsInBoundingBox(), Garfield::ComponentCST::SetRange(), SetRange(), Garfield::ComponentAnsys121::SetRangeZ(), Garfield::ComponentCST::SetRangeZ(), and UpdatePeriodicityCommon().

The documentation for this class was generated from the following files: