44 if (!this->mUseNonConstantConductivities && !this->mUseFibreOrientation)
47 c_matrix<double, SPACE_DIM, SPACE_DIM> conductivity_matrix(zero_matrix<double>(SPACE_DIM,SPACE_DIM));
49 for (
unsigned dim=0; dim<SPACE_DIM; dim++)
51 assert(this->mConstantConductivities(dim) != DBL_MAX);
52 conductivity_matrix(dim,dim) = this->mConstantConductivities(dim);
55 this->mTensors.push_back(conductivity_matrix);
59 c_matrix<double,SPACE_DIM,SPACE_DIM> orientation_matrix((identity_matrix<double>(SPACE_DIM)));
61 if (this->mUseFibreOrientation)
65 if (this->mFileReader->GetNumLinesOfData() != this->mpMesh->GetNumElements())
67 EXCEPTION(
"The size of the fibre file does not match the number of elements in the mesh");
71 if (this->mUseNonConstantConductivities)
73 if (this->mpNonConstantConductivities->size() != this->mpMesh->GetNumLocalElements())
75 EXCEPTION(
"The size of the conductivities vector does not match the number of elements in the mesh");
81 this->mTensors.reserve(this->mpMesh->GetNumLocalElements());
83 c_matrix<double, SPACE_DIM, SPACE_DIM> conductivity_matrix(zero_matrix<double>(SPACE_DIM,SPACE_DIM));
85 unsigned local_element_index = 0;
88 it != this->mpMesh->GetElementIteratorEnd();
110 if (this->mUseNonConstantConductivities)
112 for (
unsigned dim=0; dim<SPACE_DIM; dim++)
114 conductivity_matrix(dim,dim) = (*this->mpNonConstantConductivities)[local_element_index][dim];
119 for (
unsigned dim=0; dim<SPACE_DIM; dim++)
121 assert(this->mConstantConductivities(dim) != DBL_MAX);
122 conductivity_matrix(dim,dim) = this->mConstantConductivities(dim);
126 if (this->mUseFibreOrientation)
128 unsigned current_fibre_global_index = it->GetIndex();
129 this->mFileReader->GetFibreSheetAndNormalMatrix(current_fibre_global_index, orientation_matrix);
132 c_matrix<double,SPACE_DIM,SPACE_DIM> temp;
133 noalias(temp) = prod(orientation_matrix, conductivity_matrix);
134 this->mTensors.push_back( prod(temp, trans(orientation_matrix) ) );
136 local_element_index++;
138 assert(this->mTensors.size() == this->mpMesh->GetNumLocalElements());
139 assert(this->mTensors.size() == local_element_index);
141 if (this->mUseFibreOrientation)
144 this->mFileReader.reset();
148 this->mInitialised =
true;