Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Conditions.cpp
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////////////
2 //
3 // File: BoundaryConditions.cpp
4 //
5 // For more information, please see: http://www.nektar.info/
6 //
7 // The MIT License
8 //
9 // Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10 // Department of Aeronautics, Imperial College London (UK), and Scientific
11 // Computing and Imaging Institute, University of Utah (USA).
12 //
13 // License for the specific language governing rights and limitations under
14 // Permission is hereby granted, free of charge, to any person obtaining a
15 // copy of this software and associated documentation files (the "Software"),
16 // to deal in the Software without restriction, including without limitation
17 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
18 // and/or sell copies of the Software, and to permit persons to whom the
19 // Software is furnished to do so, subject to the following conditions:
20 //
21 // The above copyright notice and this permission notice shall be included
22 // in all copies or substantial portions of the Software.
23 //
24 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
25 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
26 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
27 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
28 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
29 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
30 // DEALINGS IN THE SOFTWARE.
31 //
32 // Description:
33 //
34 //
35 ////////////////////////////////////////////////////////////////////////////////
36 
39 #include <tinyxml.h>
40 
41 namespace Nektar
42 {
43  namespace SpatialDomains
44  {
46  : m_meshGraph(meshGraph),
47  m_session (pSession)
48 
49  {
50  Read(m_session->GetElement("Nektar/Conditions"));
51  }
52 
53 
55  {
56  }
57 
59  {
60  }
61 
62 
63  /**
64  *
65  */
66  void BoundaryConditions::Read(TiXmlElement *conditions)
67  {
68  ASSERTL0(conditions, "Unable to find CONDITIONS tag in file.");
69 
70  TiXmlElement *boundaryRegions = conditions->FirstChildElement("BOUNDARYREGIONS");
71 
72  if(boundaryRegions)
73  {
74  ReadBoundaryRegions(conditions);
75 
76  ReadBoundaryConditions(conditions);
77  }
78  }
79 
80 
81  /**
82  *
83  */
84  void BoundaryConditions::ReadBoundaryRegions(TiXmlElement *conditions)
85  {
86  // ensure boundary regions only read once per class definition
87  if(m_boundaryRegions.size() != 0)
88  {
89  return;
90  }
91 
92  TiXmlElement *boundaryRegions = conditions->FirstChildElement("BOUNDARYREGIONS");
93  ASSERTL0(boundaryRegions, "Unable to find BOUNDARYREGIONS block.");
94 
95  // See if we have boundary regions defined.
96  TiXmlElement *boundaryRegionsElement = boundaryRegions->FirstChildElement("B");
97 
98  while (boundaryRegionsElement)
99  {
100  /// All elements are of the form: "<B ID="#"> ... </B>", with
101  /// ? being the element type.
102  int indx;
103  int err = boundaryRegionsElement->QueryIntAttribute("ID", &indx);
104  ASSERTL0(err == TIXML_SUCCESS, "Unable to read attribute ID.");
105 
106  TiXmlNode* boundaryRegionChild = boundaryRegionsElement->FirstChild();
107  // This is primarily to skip comments that may be present.
108  // Comments appear as nodes just like elements.
109  // We are specifically looking for text in the body
110  // of the definition.
111  while(boundaryRegionChild && boundaryRegionChild->Type() != TiXmlNode::TEXT)
112  {
113  boundaryRegionChild = boundaryRegionChild->NextSibling();
114  }
115 
116  ASSERTL0(boundaryRegionChild, "Unable to read variable definition body.");
117  std::string boundaryRegionStr = boundaryRegionChild->ToText()->ValueStr();
118 
119  std::string::size_type indxBeg = boundaryRegionStr.find_first_of('[') + 1;
120  std::string::size_type indxEnd = boundaryRegionStr.find_last_of(']') - 1;
121 
122  ASSERTL0(indxBeg <= indxEnd, (std::string("Error reading boundary region definition:") + boundaryRegionStr).c_str());
123 
124  std::string indxStr = boundaryRegionStr.substr(indxBeg, indxEnd - indxBeg + 1);
125 
126  if (!indxStr.empty())
127  {
128  // Extract the composites from the string and return them in a list.
130 
131  ASSERTL0(m_boundaryRegions.count(indx) == 0,
132  "Boundary region "+indxStr+ " defined more than "
133  "once!");
134 
135  m_meshGraph->GetCompositeList(indxStr, *boundaryRegion);
136  m_boundaryRegions[indx] = boundaryRegion;
137  }
138 
139  boundaryRegionsElement = boundaryRegionsElement->NextSiblingElement("B");
140  }
141  }
142 
143 
144  /**
145  *
146  */
147  void BoundaryConditions::ReadBoundaryConditions(TiXmlElement *conditions)
148  {
149  // Protect against multiple reads.
150  if(m_boundaryConditions.size() != 0)
151  {
152  return;
153  }
154 
155  // Read REGION tags
156  TiXmlElement *boundaryConditionsElement = conditions->FirstChildElement("BOUNDARYCONDITIONS");
157  ASSERTL0(boundaryConditionsElement, "Boundary conditions must be specified.");
158 
159  TiXmlElement *regionElement = boundaryConditionsElement->FirstChildElement("REGION");
160 
161  // Read R (Robin), D (Dirichlet), N (Neumann), P (Periodic) C(Cauchy) tags
162  while (regionElement)
163  {
165 
166  int boundaryRegionID;
167  int err = regionElement->QueryIntAttribute("REF", &boundaryRegionID);
168  ASSERTL0(err == TIXML_SUCCESS, "Error reading boundary region reference.");
169 
170  // Find the boundary region corresponding to this ID.
171  std::string boundaryRegionIDStr;
172  std::ostringstream boundaryRegionIDStrm(boundaryRegionIDStr);
173  boundaryRegionIDStrm << boundaryRegionID;
174 
175  ASSERTL0(m_boundaryRegions.count(boundaryRegionID) == 1,
176  "Boundary region " + boost::lexical_cast<
177  string>(boundaryRegionID)+ " not found");
178 
179  TiXmlElement *conditionElement = regionElement->FirstChildElement();
180  std::vector<std::string> vars = m_session->GetVariables();
181 
182  while (conditionElement)
183  {
184  // Check type.
185  std::string conditionType = conditionElement->Value();
186  std::string attrData;
187 
188  // All have var specified, or else all variables are zero.
189  TiXmlAttribute *attr = conditionElement->FirstAttribute();
190 
192  std::string attrName;
193 
194  attrData = conditionElement->Attribute("VAR");
195 
196  if (!attrData.empty())
197  {
198  iter = std::find(vars.begin(), vars.end(), attrData);
199  ASSERTL0(iter != vars.end(), (std::string("Cannot find variable: ") + attrData).c_str());
200  }
201 
202  if (conditionType == "N")
203  {
204  if (attrData.empty())
205  {
206  // All variables are Neumann and are set to zero.
207  for (std::vector<std::string>::iterator varIter = vars.begin();
208  varIter != vars.end(); ++varIter)
209  {
211  (*boundaryConditions)[*varIter] = neumannCondition;
212  }
213  }
214  else
215  {
216  // Use the iterator from above, which must point to the variable.
217  attr = attr->Next();
218 
219  if (attr)
220  {
221  std::string equation, userDefined, filename;
222 
223  while(attr)
224  {
225 
226  attrName = attr->Name();
227 
228  if (attrName=="USERDEFINEDTYPE")
229  {
230 
231  // Do stuff for the user defined attribute
232  attrData = attr->Value();
233  ASSERTL0(!attrData.empty(), "USERDEFINEDTYPE attribute must have associated value.");
234 
235  // Suppose to go here?
236  m_session->SubstituteExpressions(attrData);
237 
238  userDefined = attrData;
239  }
240  else if(attrName=="VALUE")
241  {
242  ASSERTL0(attrName == "VALUE", (std::string("Unknown attribute: ") + attrName).c_str());
243 
244  attrData = attr->Value();
245  ASSERTL0(!attrData.empty(), "VALUE attribute must be specified.");
246 
247  m_session->SubstituteExpressions(attrData);
248 
249  equation = attrData;
250  }
251  else if(attrName=="FILE")
252  {
253  ASSERTL0(attrName == "FILE", (std::string("Unknown attribute: ") + attrName).c_str());
254 
255  attrData = attr->Value();
256  ASSERTL0(!attrData.empty(), "FILE attribute must be specified.");
257 
258  m_session->SubstituteExpressions(attrData);
259 
260  filename = attrData;
261  }
262  attr = attr->Next();
263  }
264  BoundaryConditionShPtr neumannCondition(MemoryManager<NeumannBoundaryCondition>::AllocateSharedPtr(m_session, equation, userDefined, filename));
265  (*boundaryConditions)[*iter] = neumannCondition;
266  }
267  else
268  {
269  // This variable's condition is zero.
271  (*boundaryConditions)[*iter] = neumannCondition;
272  }
273  }
274  }
275  else if (conditionType == "D")
276  {
277  if (attrData.empty())
278  {
279  // All variables are Dirichlet and are set to zero.
280  for (std::vector<std::string>::iterator varIter = vars.begin();
281  varIter != vars.end(); ++varIter)
282  {
284  (*boundaryConditions)[*varIter] = dirichletCondition;
285  }
286  }
287  else
288  {
289  // Use the iterator from above, which must point to the variable.
290  attr = attr->Next();
291 
292  if (attr)
293  {
294  std::string equation, userDefined, filename;
295 
296  while(attr) {
297 
298  attrName = attr->Name();
299 
300  if (attrName=="USERDEFINEDTYPE") {
301 
302  // Do stuff for the user defined attribute
303  attrData = attr->Value();
304  ASSERTL0(!attrData.empty(), "USERDEFINEDTYPE attribute must have associated value.");
305 
306  m_session->SubstituteExpressions(attrData);
307 
308  userDefined = attrData;
309  }
310  else if(attrName=="VALUE")
311  {
312  ASSERTL0(attrName == "VALUE", (std::string("Unknown attribute: ") + attrName).c_str());
313 
314  attrData = attr->Value();
315  ASSERTL0(!attrData.empty(), "VALUE attribute must have associated value.");
316 
317  m_session->SubstituteExpressions(attrData);
318 
319  equation = attrData;
320  }
321  else if(attrName=="FILE")
322  {
323  ASSERTL0(attrName == "FILE", (std::string("Unknown attribute: ") + attrName).c_str());
324 
325  attrData = attr->Value();
326  ASSERTL0(!attrData.empty(), "FILE attribute must be specified.");
327 
328  m_session->SubstituteExpressions(attrData);
329 
330  filename = attrData;
331  }
332  attr = attr->Next();
333  }
334 
335  BoundaryConditionShPtr dirichletCondition(MemoryManager<DirichletBoundaryCondition>::AllocateSharedPtr(m_session, equation, userDefined, filename));
336  (*boundaryConditions)[*iter] = dirichletCondition;
337  }
338  else
339  {
340  // This variable's condition is zero.
342  (*boundaryConditions)[*iter] = dirichletCondition;
343  }
344  }
345  }
346  else if (conditionType == "R") // Read du/dn + PRIMCOEFF u = VALUE
347  {
348  if (attrData.empty())
349  {
350  // All variables are Robin and are set to zero.
351  for (std::vector<std::string>::iterator varIter = vars.begin();
352  varIter != vars.end(); ++varIter)
353  {
355  (*boundaryConditions)[*varIter] = robinCondition;
356  }
357  }
358  else
359  {
360  // Use the iterator from above, which must
361  // point to the variable. Read the A and
362  // B attributes.
363  attr = attr->Next();
364 
365  if (attr)
366  {
367  std::string attrName1;
368  std::string attrData1;
369  std::string equation1, equation2, userDefined;
370  std::string filename;
371 
372  while(attr){
373 
374  attrName1 = attr->Name();
375 
376  if (attrName1=="USERDEFINEDTYPE") {
377 
378  // Do stuff for the user defined attribute
379  attrData1 = attr->Value();
380  ASSERTL0(!attrData1.empty(), "USERDEFINEDTYPE attribute must have associated value.");
381 
382  m_session->SubstituteExpressions(attrData1);
383 
384  userDefined = attrData1;
385 
386  }
387  else if(attrName1 == "VALUE"){
388 
389  ASSERTL0(attrName1 == "VALUE", (std::string("Unknown attribute: ") + attrName1).c_str());
390 
391  attrData1 = attr->Value();
392  ASSERTL0(!attrData1.empty(), "VALUE attributes must have associated values.");
393 
394  m_session->SubstituteExpressions(attrData1);
395 
396  equation1 = attrData1;
397 
398  attr = attr->Next();
399  ASSERTL0(attr, "Unable to read PRIMCOEFF attribute.");
400 
401  attrName1= attr->Name();
402  ASSERTL0(attrName1 == "PRIMCOEFF", (std::string("Unknown attribute: ") + attrName1).c_str());
403 
404  attrData1 = attr->Value();
405  ASSERTL0(!attrData1.empty(), "PRIMCOEFF attributes must have associated values.");
406 
407  m_session->SubstituteExpressions(attrData1);
408 
409  equation2 = attrData1;
410 
411  }
412  else if(attrName1=="FILE")
413  {
414  ASSERTL0(attrName1 == "FILE", (std::string("Unknown attribute: ") + attrName1).c_str());
415 
416  attrData1 = attr->Value();
417  ASSERTL0(!attrData1.empty(), "FILE attribute must be specified.");
418 
419  m_session->SubstituteExpressions(attrData1);
420 
421  filename = attrData1;
422  }
423  attr = attr->Next();
424 
425  }
426 
427  BoundaryConditionShPtr robinCondition(MemoryManager<RobinBoundaryCondition>::AllocateSharedPtr(m_session, equation1, equation2, userDefined, filename));
428  (*boundaryConditions)[*iter] = robinCondition;
429  }
430  else
431  {
432  // This variable's condition is zero.
434  (*boundaryConditions)[*iter] = robinCondition;
435  }
436  }
437  }
438  else if (conditionType == "P")
439  {
440  if (attrData.empty())
441  {
442  attr = attr->Next();
443 
444  if (attr)
445  {
446  attrName = attr->Name();
447 
448  ASSERTL0(attrName == "VALUE", (std::string("Unknown attribute: ") + attrName).c_str());
449 
450  attrData = attr->Value();
451  ASSERTL0(!attrData.empty(), "VALUE attribute must have associated value.");
452 
453  int beg = attrData.find_first_of("[");
454  int end = attrData.find_first_of("]");
455  std::string periodicBndRegionIndexStr = attrData.substr(beg+1,end-beg-1);
456  ASSERTL0(beg < end, (std::string("Error reading periodic boundary region definition for boundary region: ")
457  + boundaryRegionIDStrm.str()).c_str());
458 
459  vector<unsigned int> periodicBndRegionIndex;
460  bool parseGood = ParseUtils::GenerateSeqVector(periodicBndRegionIndexStr.c_str(), periodicBndRegionIndex);
461 
462  ASSERTL0(parseGood && (periodicBndRegionIndex.size()==1), (std::string("Unable to read periodic boundary condition for boundary region: ")
463  + boundaryRegionIDStrm.str()).c_str());
464 
465  BoundaryConditionShPtr periodicCondition(MemoryManager<PeriodicBoundaryCondition>::AllocateSharedPtr(periodicBndRegionIndex[0]));
466 
467  for (std::vector<std::string>::iterator varIter = vars.begin();
468  varIter != vars.end(); ++varIter)
469  {
470  (*boundaryConditions)[*varIter] = periodicCondition;
471  }
472  }
473  else
474  {
475  ASSERTL0(false, "Periodic boundary conditions should be explicitely defined");
476  }
477  }
478  else
479  {
480  // Use the iterator from above, which must point to the variable.
481  // Read the VALUE attribute. It is the next and only other attribute.
482  attr = attr->Next();
483 
484  if (attr)
485  {
486  attrName = attr->Name();
487 
488  ASSERTL0(attrName == "VALUE", (std::string("Unknown attribute: ") + attrName).c_str());
489 
490  attrData = attr->Value();
491  ASSERTL0(!attrData.empty(), "VALUE attribute must have associated value.");
492 
493  int beg = attrData.find_first_of("[");
494  int end = attrData.find_first_of("]");
495  std::string periodicBndRegionIndexStr = attrData.substr(beg+1,end-beg-1);
496  ASSERTL0(beg < end, (std::string("Error reading periodic boundary region definition for boundary region: ") + boundaryRegionIDStrm.str()).c_str());
497 
498  vector<unsigned int> periodicBndRegionIndex;
499  bool parseGood = ParseUtils::GenerateSeqVector(periodicBndRegionIndexStr.c_str(), periodicBndRegionIndex);
500 
501  ASSERTL0(parseGood && (periodicBndRegionIndex.size()==1), (std::string("Unable to read periodic boundary condition for boundary region: ") + boundaryRegionIDStrm.str()).c_str());
502 
503  BoundaryConditionShPtr periodicCondition(MemoryManager<PeriodicBoundaryCondition>::AllocateSharedPtr(periodicBndRegionIndex[0]));
504  (*boundaryConditions)[*iter] = periodicCondition;
505  }
506  else
507  {
508  ASSERTL0(false, "Periodic boundary conditions should be explicitely defined");
509  }
510  }
511  }
512  else if (conditionType == "C")
513  {
514  NEKERROR(ErrorUtil::ewarning, "Cauchy type boundary conditions not implemented.");
515  }
516 
517  conditionElement = conditionElement->NextSiblingElement();
518  }
519 
520  m_boundaryConditions[boundaryRegionID] = boundaryConditions;
521  regionElement = regionElement->NextSiblingElement("REGION");
522  }
523  }
524  }
525 }