69 EXCEPTION(
"Axisymmetric anisotropic conductivity only makes sense in 3D");
73 if (!this->mUseNonConstantConductivities && !this->mUseFibreOrientation)
76 c_matrix<double, SPACE_DIM, SPACE_DIM> conductivity_matrix(zero_matrix<double>(SPACE_DIM,SPACE_DIM));
78 for (
unsigned dim=0; dim<SPACE_DIM; dim++)
80 assert(this->mConstantConductivities(dim) != DBL_MAX);
81 conductivity_matrix(dim,dim) = this->mConstantConductivities(dim);
83 this->mTensors.push_back(conductivity_matrix);
87 c_vector<double,SPACE_DIM> fibre_vector((zero_vector<double>(SPACE_DIM)));
90 if (this->mUseFibreOrientation)
94 if (this->mFileReader->GetNumLinesOfData() != this->mpMesh->GetNumElements())
96 EXCEPTION(
"The size of the fibre file does not match the number of elements in the mesh");
100 if (this->mUseNonConstantConductivities)
102 if (this->mpNonConstantConductivities->size() != this->mpMesh->GetNumLocalElements())
104 EXCEPTION(
"The size of the conductivities vector does not match the number of elements in the mesh");
110 this->mTensors.reserve(this->mpMesh->GetNumLocalElements());
112 c_matrix<double, SPACE_DIM, SPACE_DIM> conductivity_matrix(zero_matrix<double>(SPACE_DIM,SPACE_DIM));
114 unsigned local_element_index = 0;
117 it != this->mpMesh->GetElementIteratorEnd();
145 if (this->mUseNonConstantConductivities)
147 for (
unsigned dim=0; dim<SPACE_DIM; dim++)
149 conductivity_matrix(dim,dim) = (*this->mpNonConstantConductivities)[local_element_index][dim];
154 for (
unsigned dim=0; dim<SPACE_DIM; dim++)
156 assert(this->mConstantConductivities(dim) != DBL_MAX);
157 conductivity_matrix(dim,dim) = this->mConstantConductivities(dim);
161 if (this->mUseFibreOrientation)
163 unsigned current_fibre_global_index = it->GetIndex();
164 this->mFileReader->GetFibreVector(current_fibre_global_index, fibre_vector);
167 this->mTensors.push_back( conductivity_matrix(1,1) * identity_matrix<double>(SPACE_DIM) +
168 (conductivity_matrix(0,0) - conductivity_matrix(1,1)) * outer_prod(fibre_vector,fibre_vector));
170 local_element_index++;
173 assert(this->mTensors.size() == this->mpMesh->GetNumLocalElements());
174 assert(this->mTensors.size() == local_element_index);
176 if (this->mUseFibreOrientation)
179 this->mFileReader.reset();
183 this->mInitialised =
true;