CellwiseDataGradient.cpp

00001 /*
00002 
00003 Copyright (C) University of Oxford, 2005-2011
00004 
00005 University of Oxford means the Chancellor, Masters and Scholars of the
00006 University of Oxford, having an administrative office at Wellington
00007 Square, Oxford OX1 2JD, UK.
00008 
00009 This file is part of Chaste.
00010 
00011 Chaste is free software: you can redistribute it and/or modify it
00012 under the terms of the GNU Lesser General Public License as published
00013 by the Free Software Foundation, either version 2.1 of the License, or
00014 (at your option) any later version.
00015 
00016 Chaste is distributed in the hope that it will be useful, but WITHOUT
00017 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00018 FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
00019 License for more details. The offer of Chaste under the terms of the
00020 License is subject to the License being interpreted in accordance with
00021 English Law and subject to any action against the University of Oxford
00022 being under the jurisdiction of the English Courts.
00023 
00024 You should have received a copy of the GNU Lesser General Public License
00025 along with Chaste. If not, see <http://www.gnu.org/licenses/>.
00026 
00027 */
00028 #include "CellwiseDataGradient.hpp"
00029 #include "LinearBasisFunction.hpp"
00030 
00031 template<unsigned DIM>
00032 c_vector<double, DIM>& CellwiseDataGradient<DIM>::rGetGradient(unsigned nodeIndex)
00033 {
00034     return mGradients[nodeIndex];
00035 }
00036 
00037 
00038 template<unsigned DIM>
00039 void CellwiseDataGradient<DIM>::SetupGradients()
00040 {
00041     MeshBasedCellPopulation<DIM>* p_cell_population = static_cast<MeshBasedCellPopulation<DIM>*>(&(CellwiseData<DIM>::Instance()->rGetCellPopulation()));
00042     TetrahedralMesh<DIM,DIM>& r_mesh = p_cell_population->rGetMesh();
00043 
00044     // Initialise gradients size
00045     unsigned num_nodes = p_cell_population->GetNumNodes();
00046     mGradients.resize(num_nodes, zero_vector<double>(DIM));
00047 
00048     // The constant gradients at each element
00049     std::vector<c_vector<double, DIM> > gradients_on_elements;
00050     unsigned num_elements = r_mesh.GetNumElements();
00051     gradients_on_elements.resize(num_elements, zero_vector<double>(DIM));
00052 
00053     // The number of elements containing a given node (excl ghost elements)
00054     std::vector<unsigned> num_real_elems_for_node(num_nodes, 0);
00055 
00056     for (unsigned elem_index=0; elem_index<num_elements; elem_index++)
00057     {
00058         Element<DIM,DIM>& r_elem = *(r_mesh.GetElement(elem_index));
00059 
00060         // Calculate the basis functions at any point (eg zero) in the element
00061         c_matrix<double, DIM, DIM> jacobian, inverse_jacobian;
00062         double jacobian_det;
00063         r_mesh.GetInverseJacobianForElement(elem_index, jacobian, jacobian_det, inverse_jacobian);
00064         const ChastePoint<DIM> zero_point;
00065         c_matrix<double, DIM, DIM+1> grad_phi;
00066         LinearBasisFunction<DIM>::ComputeTransformedBasisFunctionDerivatives(zero_point, inverse_jacobian, grad_phi);
00067 
00068         bool is_ghost_element = false;
00069 
00070         for (unsigned node_index=0; node_index<DIM+1; node_index++)
00071         {
00072             unsigned node_global_index = r_elem.GetNodeGlobalIndex(node_index);
00073 
00074             // Check whether ghost element
00075             if (p_cell_population->IsGhostNode(node_global_index) == true)
00076             {
00077                 is_ghost_element = true;
00078                 break;
00079             }
00080 
00081             // If no ghost element, get PDE solution
00082             CellPtr p_cell = p_cell_population->GetCellUsingLocationIndex(node_global_index);
00083             double pde_solution = CellwiseData<DIM>::Instance()->GetValue(p_cell, 0);
00084 
00085             // Interpolate gradient
00086             for (unsigned i=0; i<DIM; i++)
00087             {
00088                 gradients_on_elements[elem_index](i) += pde_solution* grad_phi(i, node_index);
00089             }
00090         }
00091 
00092         // Add gradient at element to gradient at node
00093         if (!is_ghost_element)
00094         {
00095             for (unsigned node_index=0; node_index<DIM+1; node_index++)
00096             {
00097                 unsigned node_global_index = r_elem.GetNodeGlobalIndex(node_index);
00098                 mGradients[node_global_index] += gradients_on_elements[elem_index];
00099                 num_real_elems_for_node[node_global_index]++;
00100             }
00101         }
00102     }
00103 
00104     // Divide to obtain average gradient
00105     for (typename AbstractCellPopulation<DIM>::Iterator cell_iter = p_cell_population->Begin();
00106          cell_iter != p_cell_population->End();
00107          ++cell_iter)
00108     {
00109         unsigned node_global_index = p_cell_population->GetLocationIndexUsingCell(*cell_iter);
00110 
00111         if (!num_real_elems_for_node[node_global_index] > 0)
00112         {
00113             // The node is a real node which is not in any real element
00114             // but should be connected to some cells (if more than one cell in mesh)
00115             Node<DIM>& this_node = *(p_cell_population->GetNodeCorrespondingToCell(*cell_iter));
00116 
00117             mGradients[node_global_index] = zero_vector<double>(DIM);
00118             unsigned num_real_adjacent_nodes = 0;
00119 
00120             // Get all the adjacent nodes which correspond to real cells
00121             std::set<Node<DIM>*> real_adjacent_nodes;
00122             real_adjacent_nodes.clear();
00123 
00124             // First loop over containing elements
00125             for (typename Node<DIM>::ContainingElementIterator element_iter = this_node.ContainingElementsBegin();
00126                  element_iter != this_node.ContainingElementsEnd();
00127                  ++element_iter)
00128             {
00129                 // Then loop over nodes therein
00130                 Element<DIM,DIM>& r_adjacent_elem = *(r_mesh.GetElement(*element_iter));
00131                 for (unsigned local_node_index=0; local_node_index<DIM+1; local_node_index++)
00132                 {
00133                     unsigned adjacent_node_global_index = r_adjacent_elem.GetNodeGlobalIndex(local_node_index);
00134 
00135                     // If not a ghost node and not the node we started with
00136                     if ( p_cell_population->IsGhostNode(adjacent_node_global_index)==false && adjacent_node_global_index != node_global_index )
00137                     {
00138 
00139                         // Calculate the contribution of gradient from this node
00140                         Node<DIM>& adjacent_node = *(r_mesh.GetNode(adjacent_node_global_index));
00141 
00142                         double this_cell_concentration = CellwiseData<DIM>::Instance()->GetValue(*cell_iter, 0);
00143                         CellPtr p_adjacent_cell = p_cell_population->GetCellUsingLocationIndex(adjacent_node_global_index);
00144                         double adjacent_cell_concentration = CellwiseData<DIM>::Instance()->GetValue(p_adjacent_cell, 0);
00145 
00146                         c_vector<double, DIM> gradient_contribution = zero_vector<double>(DIM);
00147 
00148                         if (fabs(this_cell_concentration-adjacent_cell_concentration) > 100*DBL_EPSILON)
00149                         {
00150                             c_vector<double, DIM> edge_vector = r_mesh.GetVectorFromAtoB(this_node.rGetLocation(), adjacent_node.rGetLocation());
00151                             double norm_edge_vector = norm_2(edge_vector);
00152                             gradient_contribution = edge_vector
00153                                                         * (adjacent_cell_concentration - this_cell_concentration)
00154                                                         / (norm_edge_vector * norm_edge_vector);
00155                         }
00156 
00157                         mGradients[node_global_index] += gradient_contribution;
00158                         num_real_adjacent_nodes++;
00159                     }
00160                 }
00161             }
00162             mGradients[node_global_index] /= num_real_adjacent_nodes;
00163         }
00164         else
00165         {
00166             mGradients[node_global_index] /= num_real_elems_for_node[node_global_index];
00167         }
00168     }
00169 }
00170 
00171 
00173 // Explicit instantiation
00175 
00176 template class CellwiseDataGradient<1>;
00177 template class CellwiseDataGradient<2>;
00178 template class CellwiseDataGradient<3>;

Generated on Tue May 31 14:31:40 2011 for Chaste by  doxygen 1.5.5