Chaste  Release::2018.1
ContinuumMechanicsProblemDefinition.cpp
1 /*
2 
3 Copyright (c) 2005-2018, University of Oxford.
4 All rights reserved.
5 
6 University of Oxford means the Chancellor, Masters and Scholars of the
7 University of Oxford, having an administrative office at Wellington
8 Square, Oxford OX1 2JD, UK.
9 
10 This file is part of Chaste.
11 
12 Redistribution and use in source and binary forms, with or without
13 modification, are permitted provided that the following conditions are met:
14  * Redistributions of source code must retain the above copyright notice,
15  this list of conditions and the following disclaimer.
16  * Redistributions in binary form must reproduce the above copyright notice,
17  this list of conditions and the following disclaimer in the documentation
18  and/or other materials provided with the distribution.
19  * Neither the name of the University of Oxford nor the names of its
20  contributors may be used to endorse or promote products derived from this
21  software without specific prior written permission.
22 
23 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
24 AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
25 IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
26 ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
27 LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
28 CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
29 GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
30 HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
31 LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
32 OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
33 
34 */
35 
36 #include <limits>
37 
38 #include "ContinuumMechanicsProblemDefinition.hpp"
39 #include "AbstractIncompressibleMaterialLaw.hpp"
40 #include "AbstractCompressibleMaterialLaw.hpp"
41 
42 
43 template<unsigned DIM>
44 const double ContinuumMechanicsProblemDefinition<DIM>::FREE = std::numeric_limits<double>::max();
45 
46 template<unsigned DIM>
48  : mrMesh(rMesh),
49  mDensity(1.0),
50  mBodyForceType(CONSTANT_BODY_FORCE),
51  mConstantBodyForce(zero_vector<double>(DIM)),
52  mTractionBoundaryConditionType(NO_TRACTIONS),
53  mVerboseDuringSolve(false)
54 {
55 }
56 
57 template<unsigned DIM>
59 {
60  assert(density>0.0);
61  mDensity = density;
62 }
63 
64 template<unsigned DIM>
66 {
67  return mDensity;
68 }
69 
70 template<unsigned DIM>
71 void ContinuumMechanicsProblemDefinition<DIM>::SetBodyForce(c_vector<double,DIM> bodyForce)
72 {
73  mBodyForceType = CONSTANT_BODY_FORCE;
74  mConstantBodyForce = bodyForce;
75 }
76 
77 template<unsigned DIM>
78 void ContinuumMechanicsProblemDefinition<DIM>::SetBodyForce(c_vector<double,DIM> (*pFunction)(c_vector<double,DIM>& rX, double t))
79 {
80  mBodyForceType = FUNCTIONAL_BODY_FORCE;
81  mpBodyForceFunction = pFunction;
82 }
83 
84 template<unsigned DIM>
86 {
87  return mBodyForceType;
88 }
89 
90 template<unsigned DIM>
92 {
93  assert(mBodyForceType==CONSTANT_BODY_FORCE);
94  return mConstantBodyForce;
95 }
96 
97 template<unsigned DIM>
98 c_vector<double,DIM> ContinuumMechanicsProblemDefinition<DIM>::EvaluateBodyForceFunction(c_vector<double,DIM>& rX, double t)
99 {
100  assert(mBodyForceType==FUNCTIONAL_BODY_FORCE);
101  return (*mpBodyForceFunction)(rX,t);
102 }
103 
104 template<unsigned DIM>
105 c_vector<double,DIM> ContinuumMechanicsProblemDefinition<DIM>::GetBodyForce(c_vector<double,DIM>& rX, double t)
106 {
107  switch(mBodyForceType)
108  {
109  case CONSTANT_BODY_FORCE:
110  {
111  return mConstantBodyForce;
112  }
113  case FUNCTIONAL_BODY_FORCE:
114  {
115  return (*mpBodyForceFunction)(rX,t);
116  }
117  default:
119  }
120 }
121 
122 template<unsigned DIM>
124 {
125  return mTractionBoundaryConditionType;
126 }
127 
128 template<unsigned DIM>
130  std::vector<c_vector<double,DIM> >& rElementwiseTractions)
131 {
132 
133  assert(rTractionBoundaryElements.size()==rElementwiseTractions.size());
134  mTractionBoundaryConditionType = ELEMENTWISE_TRACTION;
135  mTractionBoundaryElements = rTractionBoundaryElements;
136  mElementwiseTractions = rElementwiseTractions;
137 }
138 
139 template<unsigned DIM>
141  c_vector<double,DIM> (*pFunction)(c_vector<double,DIM>& rX, double t))
142 {
143  mTractionBoundaryConditionType=FUNCTIONAL_TRACTION;
144  mTractionBoundaryElements = rTractionBoundaryElements;
145  mpTractionBoundaryConditionFunction = pFunction;
146 }
147 
148 template<unsigned DIM>
150  double normalPressure)
151 {
152  mTractionBoundaryConditionType = PRESSURE_ON_DEFORMED;
153  mTractionBoundaryElements = rTractionBoundaryElements;
154  mNormalPressure = normalPressure;
155  mOriginalNormalPressure = normalPressure;
156 
157 }
158 
159 template<unsigned DIM>
161  double (*pFunction)(double t))
162 {
163  mTractionBoundaryConditionType = FUNCTIONAL_PRESSURE_ON_DEFORMED;
164  mTractionBoundaryElements = rTractionBoundaryElements;
165  mpNormalPressureFunction = pFunction;
166 }
167 
168 template<unsigned DIM>
169 void ContinuumMechanicsProblemDefinition<DIM>::SetZeroDirichletNodes(std::vector<unsigned>& rZeroDirichletNodes)
170 {
171  mDirichletNodes = rZeroDirichletNodes;
172 
173  for (unsigned i=0; i<mDirichletNodes.size(); i++)
174  {
175  assert(mDirichletNodes[i] < mrMesh.GetNumNodes());
176  }
177 
178  mDirichletNodeValues.clear();
179  for (unsigned i=0; i<mDirichletNodes.size(); i++)
180  {
181  mDirichletNodeValues.push_back(zero_vector<double>(DIM));
182  }
183 }
184 
185 template<unsigned DIM>
187 {
188  return mDirichletNodes;
189 }
190 
191 template<unsigned DIM>
193 {
194  return mDirichletNodeValues;
195 }
196 
197 template<unsigned DIM>
199 {
200  return mTractionBoundaryElements;
201 }
202 
203 template<unsigned DIM>
205 {
206  assert(mTractionBoundaryConditionType==ELEMENTWISE_TRACTION);
207  return mElementwiseTractions;
208 }
209 
210 template<unsigned DIM>
212 {
213  assert(mTractionBoundaryConditionType==PRESSURE_ON_DEFORMED);
214  return mNormalPressure;
215 }
216 
217 template<unsigned DIM>
218 c_vector<double,DIM> ContinuumMechanicsProblemDefinition<DIM>::EvaluateTractionFunction(c_vector<double,DIM>& rX, double t)
219 {
220  assert(mTractionBoundaryConditionType==FUNCTIONAL_TRACTION);
221  return (*mpTractionBoundaryConditionFunction)(rX,t);
222 }
223 
224 template<unsigned DIM>
226 {
227  assert(mTractionBoundaryConditionType==FUNCTIONAL_PRESSURE_ON_DEFORMED);
228  return (*mpNormalPressureFunction)(t);
229 }
230 
231 template<unsigned DIM>
233 {
234  assert(mTractionBoundaryConditionType==PRESSURE_ON_DEFORMED);
235  mNormalPressure = mOriginalNormalPressure*scaleFactor;
236 }
237 
238 template<unsigned DIM>
240 {
241  if (mDirichletNodes.size() == 0)
242  {
243  EXCEPTION("No Dirichlet boundary conditions (eg fixed displacement or fixed flow) have been set");
244  }
245 }
246 
247 // Explicit instantiation
std::vector< BoundaryElement< DIM-1, DIM > * > & rGetTractionBoundaryElements()
ContinuumMechanicsProblemDefinition(AbstractTetrahedralMesh< DIM, DIM > &rMesh)
void SetTractionBoundaryConditions(std::vector< BoundaryElement< DIM-1, DIM > * > &rTractionBoundaryElements, std::vector< c_vector< double, DIM > > &rElementwiseTractions)
c_vector< double, DIM > EvaluateTractionFunction(c_vector< double, DIM > &rX, double t)
#define EXCEPTION(message)
Definition: Exception.hpp:143
TractionBoundaryConditionType GetTractionBoundaryConditionType()
void SetZeroDirichletNodes(std::vector< unsigned > &rZeroDirichletNodes)
std::vector< c_vector< double, DIM > > & rGetElementwiseTractions()
#define NEVER_REACHED
Definition: Exception.hpp:206
c_vector< double, DIM > GetBodyForce(c_vector< double, DIM > &rX, double t=0.0)
void SetBodyForce(c_vector< double, DIM > bodyForce)
std::vector< c_vector< double, DIM > > & rGetDirichletNodeValues()
void SetApplyNormalPressureOnDeformedSurface(std::vector< BoundaryElement< DIM-1, DIM > * > &rTractionBoundaryElements, double normalPressure)
c_vector< double, DIM > EvaluateBodyForceFunction(c_vector< double, DIM > &rX, double t)