1 /*=========================================================================
2 *
3 * Copyright Insight Software Consortium
4 *
5 * Licensed under the Apache License, Version 2.0 (the "License");
6 * you may not use this file except in compliance with the License.
7 * You may obtain a copy of the License at
8 *
9 * http://www.apache.org/licenses/LICENSE-2.0.txt
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 *
17 *=========================================================================*/
18
19 #include "metaFEMObject.h"
20
21 #include <stdio.h>
22 #include <ctype.h>
23 #include <string>
24
25 #if (METAIO_USE_NAMESPACE)
26 namespace METAIO_NAMESPACE {
27 #endif
28
29 FEMObjectNode::
FEMObjectNode(int dim)30 FEMObjectNode(int dim)
31 {
32 m_Dim = dim;
33 m_GN = -1;
34 m_X = new float[m_Dim];
35 for(unsigned int i=0;i<m_Dim;i++)
36 {
37 m_X[i] = 0;
38 }
39 }
40
41 FEMObjectNode::
~FEMObjectNode()42 ~FEMObjectNode()
43 {
44 delete []m_X;
45 }
46
47 FEMObjectElement::
FEMObjectElement(int dim)48 FEMObjectElement(int dim)
49 {
50 m_Dim = dim;
51 m_GN = -1;
52 m_NodesId = new int[m_Dim];
53 for(unsigned int i=0;i<m_Dim;i++)
54 {
55 m_NodesId[i] = -1;
56 }
57 }
58
59 FEMObjectElement::
~FEMObjectElement()60 ~FEMObjectElement()
61 {
62 delete []m_NodesId;
63 }
64
FEMObjectLoad()65 FEMObjectLoad::FEMObjectLoad()
66 {
67 }
68
~FEMObjectLoad()69 FEMObjectLoad::~FEMObjectLoad()
70 {
71 for(METAIO_STL::vector<FEMObjectMFCTerm *>::iterator it = this->m_LHS.begin();
72 it != this->m_LHS.end();
73 ++it)
74 {
75 delete (*it);
76 }
77 this->m_LHS.clear();
78 this->m_RHS.clear();
79 this->m_ForceMatrix.clear();
80 this->m_ForceVector.clear();
81 }
82
83 //
84 // MetaFEMObject Constructors
85 //
86 MetaFEMObject::
MetaFEMObject()87 MetaFEMObject()
88 :MetaObject()
89 {
90 if(META_DEBUG) METAIO_STREAM::cout << "MetaFEMObject()" << METAIO_STREAM::endl;
91
92 Clear();
93
94 this->m_ClassNameList.push_back("Node");
95 this->m_ClassNameList.push_back("MaterialLinearElasticity");
96 this->m_ClassNameList.push_back("Element2DC0LinearLineStress");
97 this->m_ClassNameList.push_back("Element2DC1Beam");
98 this->m_ClassNameList.push_back("Element2DC0LinearTriangularMembrane");
99 this->m_ClassNameList.push_back("Element2DC0LinearTriangularStrain");
100 this->m_ClassNameList.push_back("Element2DC0LinearTriangularStress");
101 this->m_ClassNameList.push_back("Element2DC0LinearQuadrilateralMembrane");
102 this->m_ClassNameList.push_back("Element2DC0LinearQuadrilateralStrain");
103 this->m_ClassNameList.push_back("Element2DC0LinearQuadrilateralStress");
104 this->m_ClassNameList.push_back("Element2DC0QuadraticTriangularStress");
105 this->m_ClassNameList.push_back("Element2DC0QuadraticTriangularStrain");
106 this->m_ClassNameList.push_back("Element3DC0LinearHexahedronMembrane");
107 this->m_ClassNameList.push_back("Element3DC0LinearHexahedronStrain");
108 this->m_ClassNameList.push_back("Element3DC0LinearTetrahedronMembrane");
109 this->m_ClassNameList.push_back("Element3DC0LinearTetrahedronStrain");
110 this->m_ClassNameList.push_back("LoadBC");
111 this->m_ClassNameList.push_back("LoadBCMFC");
112 this->m_ClassNameList.push_back("LoadNode");
113 this->m_ClassNameList.push_back("LoadEdge");
114 this->m_ClassNameList.push_back("LoadGravConst");
115 this->m_ClassNameList.push_back("LoadLandmark");
116 this->m_ClassNameList.push_back("LoadPoint");
117 this->m_ElementDataFileName = "LOCAL";
118 }
119
120 //
121 MetaFEMObject::
MetaFEMObject(const char * _headerName)122 MetaFEMObject(const char *_headerName)
123 :MetaObject()
124 {
125 if(META_DEBUG)
126 {
127 METAIO_STREAM::cout << "MetaFEMObject()" << METAIO_STREAM::endl;
128 }
129 Clear();
130 Read(_headerName);
131 this->m_ElementDataFileName = "LOCAL";
132 }
133
134 //
135 MetaFEMObject::
MetaFEMObject(const MetaFEMObject * _mesh)136 MetaFEMObject(const MetaFEMObject *_mesh)
137 : MetaObject()
138 {
139 if(META_DEBUG)
140 {
141 METAIO_STREAM::cout << "MetaFEMObject()" << METAIO_STREAM::endl;
142 }
143 Clear();
144 CopyInfo(_mesh);
145 }
146
147
148
149 //
150 MetaFEMObject::
MetaFEMObject(unsigned int dim)151 MetaFEMObject(unsigned int dim)
152 :MetaObject(dim)
153 {
154 if(META_DEBUG)
155 {
156 METAIO_STREAM::cout << "MetaFEMObject()" << METAIO_STREAM::endl;
157 }
158 Clear();
159 this->m_ElementDataFileName = "LOCAL";
160 }
161
162 /** Destructor */
163 MetaFEMObject::
~MetaFEMObject()164 ~MetaFEMObject()
165 {
166 // Delete the list of pointers to Nodes.
167 NodeListType::iterator it_Node = m_NodeList.begin();
168 while(it_Node != m_NodeList.end())
169 {
170 FEMObjectNode* Node = *it_Node;
171 ++it_Node;
172 delete Node;
173 }
174 // Delete the list of pointers to Materials.
175 MaterialListType::iterator it_Material = m_MaterialList.begin();
176 while(it_Material != m_MaterialList.end())
177 {
178 FEMObjectMaterial* Material = *it_Material;
179 ++it_Material;
180 delete Material;
181 }
182
183 // Delete the list of pointers to Elements.
184 ElementListType::iterator it_Element = m_ElementList.begin();
185 while(it_Element != m_ElementList.end())
186 {
187 FEMObjectElement* Element = *it_Element;
188 ++it_Element;
189 delete Element;
190 }
191
192 // Delete the list of pointers to Loads.
193 LoadListType::iterator it_Load = m_LoadList.begin();
194 while(it_Load != m_LoadList.end())
195 {
196 FEMObjectLoad* Load = *it_Load;
197 ++it_Load;
198 delete Load;
199 }
200
201 M_Destroy();
202 }
203
204 //
205 void MetaFEMObject::
PrintInfo() const206 PrintInfo() const
207 {
208 MetaObject::PrintInfo();
209 }
210
211 void MetaFEMObject::
CopyInfo(const MetaObject * _object)212 CopyInfo(const MetaObject * _object)
213 {
214 MetaObject::CopyInfo(_object);
215 }
216
217 /** Clear FEMObject information */
218 void MetaFEMObject::
Clear(void)219 Clear(void)
220 {
221 if(META_DEBUG)
222 {
223 METAIO_STREAM::cout << "MetaFEMObject: Clear" << METAIO_STREAM::endl;
224 }
225 MetaObject::Clear();
226 if(META_DEBUG)
227 {
228 METAIO_STREAM::cout << "MetaFEMObject: Clear: m_NPoints" << METAIO_STREAM::endl;
229 }
230
231 // Delete the list of pointers to Nodes.
232 NodeListType::iterator it_Node = m_NodeList.begin();
233 while(it_Node != m_NodeList.end())
234 {
235 FEMObjectNode* Node = *it_Node;
236 ++it_Node;
237 delete Node;
238 }
239
240 // Delete the list of pointers to Elements.
241 ElementListType::iterator it_Element = m_ElementList.begin();
242 while(it_Element != m_ElementList.end())
243 {
244 FEMObjectElement* Element = *it_Element;
245 ++it_Element;
246 delete Element;
247 }
248
249 // Delete the list of pointers to Loads.
250 LoadListType::iterator it_Load = m_LoadList.begin();
251 while(it_Load != m_LoadList.end())
252 {
253 FEMObjectLoad* Load = *it_Load;
254 ++it_Load;
255 delete Load;
256 }
257
258 // Delete the list of pointers to Materials.
259 MaterialListType::iterator it_Material = m_MaterialList.begin();
260 while(it_Material != m_MaterialList.end())
261 {
262 FEMObjectMaterial* Material = *it_Material;
263 ++it_Material;
264 delete Material;
265 }
266
267 m_NodeList.clear();
268 m_ElementList.clear();
269 m_LoadList.clear();
270 m_MaterialList.clear();
271 }
272
273 /** Destroy tube information */
274 void MetaFEMObject::
M_Destroy(void)275 M_Destroy(void)
276 {
277 MetaObject::M_Destroy();
278 }
279
280 /** Set Read fields */
281 void MetaFEMObject::
M_SetupReadFields(void)282 M_SetupReadFields(void)
283 {
284 if(META_DEBUG)
285 {
286 METAIO_STREAM::cout << "MetaFEMObject: M_SetupReadFields" << METAIO_STREAM::endl;
287 }
288
289 MetaObject::M_SetupReadFields();
290
291 MET_FieldRecordType * mF;
292
293 mF = new MET_FieldRecordType;
294 MET_InitWriteField(mF, "ElementDataFile", MET_STRING, true);
295 mF->required = true;
296 mF->terminateRead = true;
297 m_Fields.push_back(mF);
298
299 }
300
301 void MetaFEMObject::
M_SetupWriteFields(void)302 M_SetupWriteFields(void)
303 {
304 strcpy(m_ObjectTypeName,"FEMObject");
305 MetaObject::M_SetupWriteFields();
306
307 MET_FieldRecordType * mF;
308
309 mF = new MET_FieldRecordType;
310 MET_InitWriteField(mF, "ElementDataFile", MET_STRING,
311 m_ElementDataFileName.length(),
312 m_ElementDataFileName.c_str());
313 mF->terminateRead = true;
314 m_Fields.push_back(mF);
315 }
316
317
318
319 bool MetaFEMObject::
M_Read(void)320 M_Read(void)
321 {
322 if(META_DEBUG)
323 {
324 METAIO_STREAM::cout << "MetaFEMObject: M_Read: Loading Header" << METAIO_STREAM::endl;
325 }
326
327 if(!MetaObject::M_Read())
328 {
329 METAIO_STREAM::cout << "MetaFEMObject: M_Read: Error parsing file" << METAIO_STREAM::endl;
330 return false;
331 }
332
333 if(META_DEBUG)
334 {
335 METAIO_STREAM::cout << "MetaFEMObject: M_Read: Parsing Header" << METAIO_STREAM::endl;
336 }
337
338 // currently reader handles only ASCII data
339 if(m_BinaryData)
340 {
341 METAIO_STREAM::cout << "MetaFEMObject: M_Read: Data content should be in ASCII format" << METAIO_STREAM::endl;
342 return false;
343 }
344
345 // we read 1)node, 2) material, 3) element and 4) load and the input should be in the afore
346 // mentioned order.
347
348 int segment_read = 0; // to keep track of what is being read. corresponds to enumerated types in the header file
349 /* then we start reading objects from stream */
350 do
351 {
352 // local variables
353 std::streampos l(0);
354 char buf[256];
355 std::string s;
356 std::string::size_type b, e;
357 bool clID;
358 std::string errorMessage;
359
360 if(segment_read > 3)
361 {
362 this->SkipWhiteSpace();
363 return true; // end of FEM segment in spatial object reader.
364 }
365
366 l = this->m_ReadStream->tellg(); // remember the stream position
367 this->SkipWhiteSpace(); // skip comments and whitespaces
368 if ( this->m_ReadStream->eof() )
369 {
370 return 0; // end of stream. all was good
371 }
372 char c;
373 if ( ( c = static_cast<char>(this->m_ReadStream->get()) ) != '<' )
374 {
375 std::string rest;
376 std::getline(*this->m_ReadStream, rest);
377 errorMessage = "Expected < token not found. Instead found '";
378 errorMessage += c;
379 errorMessage += "'.\nRest of line is '";
380 errorMessage += rest;
381 errorMessage += "'.\n";
382 METAIO_STREAM::cout << errorMessage << METAIO_STREAM::endl;
383 return false; // the file is not in proper format
384 }
385 this->m_ReadStream->getline(buf, 256, '>'); // read up to 256 characters until '>' is reached.
386 // we read and discard the '>'
387 s = std::string(buf);
388
389 // get rid of the whitespaces in front of and the back of token
390 b = s.find_first_not_of(this->whitespaces); // end of
391 // whitespaces
392 // in the
393 // beginning
394 if ( ( e = s.find_first_of(this->whitespaces, b) ) == std::string::npos ) //
395 // beginning
396 // of
397 // whitespaces
398 // at the
399 // end
400 {
401 e = s.size();
402 }
403 s = s.substr(b, e - b);
404
405 if ( s == "END" )
406 {
407 /*
408 * We can ignore this token. Start again by reading the next object.
409 */
410 segment_read++;
411 }
412 else
413 {
414 clID = this->IsClassNamePresent(s); // obtain the class ID from FEMObjectFactory
415 if ( !clID)
416 {
417 errorMessage = s;
418 errorMessage += " is not a valid FEM data type";
419 errorMessage += "'.";
420 METAIO_STREAM::cout << errorMessage << METAIO_STREAM::endl;
421 return false; // class not found
422 }
423 /*
424 * Now we have to read additional data, which is
425 * specific to the class of object we just created
426 */
427
428 switch ( segment_read )
429 {
430 case NODE :
431 this->M_Read_Node();
432 break;
433 case MATERIAL :
434 this->M_Read_Material(s);
435 break;
436 case ELEMENT :
437 this->M_Read_Element(s);
438 break;
439 case LOAD :
440 this->M_Read_Load(s);
441 break;
442 }
443 }
444 }while(segment_read <= 3); // end of FEM segment in spatial object reader.
445 return true;
446 }
447
448 bool MetaFEMObject::
M_Write(void)449 M_Write(void)
450 {
451 if(!MetaObject::M_Write())
452 {
453 METAIO_STREAM::cout << "MetaFEMObject: M_Write: Error parsing file" << METAIO_STREAM::endl;
454 return false;
455 }
456
457 NodeListType::iterator it_Node = m_NodeList.begin();
458 while(it_Node != m_NodeList.end())
459 {
460 FEMObjectNode* Node = *it_Node;
461 this->M_Write_Node(Node);
462 ++it_Node;
463 }
464 *this->m_WriteStream << "\n<END> % End of nodes\n\n";
465
466 MaterialListType::iterator it_Material = m_MaterialList.begin();
467 while(it_Material != m_MaterialList.end())
468 {
469 FEMObjectMaterial* Material = *it_Material;
470 this->M_Write_Material(Material);
471 ++it_Material;
472 }
473 *this->m_WriteStream << "\n<END> % End of material definition\n\n";
474
475
476 ElementListType::iterator it_Element = m_ElementList.begin();
477 while(it_Element != m_ElementList.end())
478 {
479 FEMObjectElement* Element = *it_Element;
480 this->M_Write_Element(Element);
481 ++it_Element;
482 }
483 *this->m_WriteStream << "\n<END> % End of element definition\n\n";
484
485 LoadListType::iterator it_Load = m_LoadList.begin();
486 while(it_Load != m_LoadList.end())
487 {
488 FEMObjectLoad* Load = *it_Load;
489 this->M_Write_Load(Load);
490 ++it_Load;
491 }
492 *this->m_WriteStream << "\n<END> % End of load definition\n\n";
493
494 return true;
495 }
496
M_Write_Node(FEMObjectNode * Node)497 void MetaFEMObject::M_Write_Node(FEMObjectNode *Node)
498 {
499 // first write the class name
500 *this->m_WriteStream << '<' << "Node" << ">\n";
501
502 // then the global object number
503 *this->m_WriteStream << "\t" << Node->m_GN << "\t% Global object number\n";
504
505 /* write co-ordinate values */
506 *this->m_WriteStream << "\t" << Node->m_Dim;
507 for (unsigned int i=0; i<Node->m_Dim; i++)
508 {
509 *this->m_WriteStream << " " << Node->m_X[i];
510 }
511 *this->m_WriteStream << "\t% Node coordinates" << "\n";
512 }
513
M_Write_Material(FEMObjectMaterial * Material)514 void MetaFEMObject::M_Write_Material(FEMObjectMaterial *Material)
515 {
516 if(std::string(Material->m_MaterialName) == "MaterialLinearElasticity")
517 {
518 *this->m_WriteStream << '<' << "MaterialLinearElasticity" << ">\n";
519 *this->m_WriteStream << "\t" << Material->m_GN << "\t% Global object number\n";
520 *this->m_WriteStream << "\tE : " << Material->E << "\t% Young modulus\n";
521 *this->m_WriteStream << "\tA : " << Material->A << "\t% Beam crossection area\n";
522 *this->m_WriteStream << "\tI : " << Material->I << "\t% Moment of inertia\n";
523 *this->m_WriteStream << "\tnu : " << Material->nu << "\t% Poisson's ratio\n";
524 *this->m_WriteStream << "\th : " << Material->h << "\t% Plate thickness\n";
525 *this->m_WriteStream << "\tRhoC : " << Material->RhoC << "\t% Density times capacity\n";
526 *this->m_WriteStream << "\tEND:\t% End of material definition\n";
527 }
528 }
529
M_Write_Element(FEMObjectElement * Element)530 void MetaFEMObject::M_Write_Element(FEMObjectElement *Element)
531 {
532 *this->m_WriteStream << '<' << Element->m_ElementName << ">\n";
533 *this->m_WriteStream << "\t" << Element->m_GN << "\t% Global object number\n";
534 unsigned int numNodes = Element->m_NumNodes;
535 for ( unsigned int p = 0; p < numNodes; p++ )
536 {
537 *this->m_WriteStream << "\t" << Element->m_NodesId[p] << "\t% Node #" << ( p + 1 ) << " ID\n";
538 }
539 *this->m_WriteStream << "\t" << Element->m_MaterialGN << "\t% Material ID\n";
540 }
541
M_Write_Load(FEMObjectLoad * Load)542 void MetaFEMObject::M_Write_Load(FEMObjectLoad *Load)
543 {
544 *this->m_WriteStream << '<' << Load->m_LoadName << ">\n";
545 *this->m_WriteStream << "\t" << Load->m_GN << "\t% Global object number\n";
546
547 // write according to the load type
548 if(std::string(Load->m_LoadName) == "LoadBC")
549 {
550 *this->m_WriteStream << "\t" << Load->m_ElementGN << "\t% GN of element" << "\n";
551 *this->m_WriteStream << "\t" << Load->m_DOF << "\t% DOF# in element" << "\n";
552 *this->m_WriteStream << "\t" << Load->m_NumRHS;
553 for (int i=0; i<Load->m_NumRHS; i++)
554 {
555 *this->m_WriteStream << " " << Load->m_RHS[i];
556 }
557 *this->m_WriteStream << "\t% value of the fixed DOF" << "\n";
558 return;
559 }
560
561 if(std::string(Load->m_LoadName) == "LoadNode")
562 {
563 *this->m_WriteStream << "\t" << Load->m_ElementGN << "\t% GN of element" << "\n";
564 *this->m_WriteStream << "\t" << Load->m_NodeNumber << " " << "\t% Point number within the element\n";
565 *this->m_WriteStream << "\t" << Load->m_Dim;
566 for (int i=0; i<Load->m_Dim; i++)
567 {
568 *this->m_WriteStream << " " << Load->m_ForceVector[i];
569 }
570 *this->m_WriteStream << "\t% Force vector (first number is the size of a vector)\n";
571 return;
572 }
573
574 if(std::string(Load->m_LoadName) == "LoadBCMFC")
575 {
576 /** write the number of DOFs affected by this MFC */
577 *this->m_WriteStream << "\t" << Load->m_NumLHS << "\t% Number of DOFs in this MFC" << std::endl;
578
579 /** write each term */
580 *this->m_WriteStream << "\t %==>\n";
581 for ( int i=0; i<Load->m_NumLHS; i++ )
582 {
583 FEMObjectMFCTerm *mfcTerm = dynamic_cast< FEMObjectMFCTerm * > (&*Load->m_LHS[i]);
584 *this->m_WriteStream << "\t " <<mfcTerm->m_ElementGN << "\t% GN of element" << std::endl;
585 *this->m_WriteStream << "\t " << mfcTerm->m_DOF << "\t% DOF# in element" << std::endl;
586 *this->m_WriteStream << "\t " << mfcTerm->m_Value << "\t% weight" << std::endl;
587 *this->m_WriteStream << "\t %==>\n";
588 }
589
590 /** write the rhs */
591 *this->m_WriteStream << "\t" << Load->m_NumRHS;
592 for ( int i=0; i<Load->m_NumRHS; i++ )
593 {
594 *this->m_WriteStream << " " << Load->m_RHS[i];
595 }
596 *this->m_WriteStream << "\t% rhs of MFC" << std::endl;
597 return;
598 }
599
600 if(std::string(Load->m_LoadName) == "LoadEdge")
601 {
602 *this->m_WriteStream << "\t" << Load->m_ElementGN << "\t% GN of the element on which the load acts" << "\n";
603 /** ... edge number */
604 *this->m_WriteStream << "\t" << Load->m_EdgeNumber << "\t% Edge number" << "\n";
605
606 /** ... force matrix */
607 size_t numRows = Load->m_ForceMatrix.size();
608 size_t numCols = Load->m_ForceMatrix[0].size();
609
610 *this->m_WriteStream << "\t" << numRows << "\t% # rows in force matrix" << "\n";
611 *this->m_WriteStream << "\t" << numCols << "\t% # cols in force matrix" << "\n";
612 *this->m_WriteStream << "\t% force matrix\n";
613 for ( size_t i = 0; i < numRows; i++ )
614 {
615 *this->m_WriteStream << "\t";
616 METAIO_STL::vector<float> F = Load->m_ForceMatrix[i];
617 for ( size_t j = 0; j < numCols; j++ )
618 {
619 *this->m_WriteStream << F[j] << " ";
620 }
621 *this->m_WriteStream << "\n";
622 }
623 return;
624 }
625
626 if(std::string(Load->m_LoadName) == "LoadGravConst")
627 {
628 /** Write the list of element global numbers */
629 if ( Load->m_NumElements > 0 )
630 {
631 *this->m_WriteStream << "\t" << Load->m_NumElements;
632 *this->m_WriteStream << "\t% # of elements on which the load acts" << std::endl;
633 *this->m_WriteStream << "\t";
634 for ( int i = 0; i < Load->m_NumElements; i++ )
635 {
636 *this->m_WriteStream << Load->m_Elements[i] << " ";
637 }
638 *this->m_WriteStream << "\t% GNs of elements" << std::endl;
639 }
640 else
641 {
642 *this->m_WriteStream << "\t-1\t% Load acts on all elements" << std::endl;
643 }
644 /** then write the actual data force vector */
645 *this->m_WriteStream << "\t" << Load->m_Dim << "\t% Size of the gravity force vector\n";
646 for (int i=0; i<Load->m_Dim; i++)
647 {
648 *this->m_WriteStream << "\t" << Load->m_ForceVector[i] ;
649 }
650 *this->m_WriteStream << "\t% Gravity force vector\n";
651 }
652
653 if(std::string(Load->m_LoadName) == "LoadLandmark")
654 {
655 // print undeformed coordinates
656 size_t dim = Load->m_Undeformed.size();
657
658 *this->m_WriteStream << "\t" << dim;
659 for ( size_t i = 0; i < dim; i++ )
660 {
661 *this->m_WriteStream << Load->m_Undeformed[i] << " ";
662 }
663 *this->m_WriteStream << "\t % Dimension , undeformed state local coordinates";
664 *this->m_WriteStream << "\n";
665
666 // print deformed coordinates
667 *this->m_WriteStream << "\t" << dim;
668 for ( size_t i = 0; i < dim; i++ )
669 {
670 *this->m_WriteStream << Load->m_Deformed[i] << " ";
671 }
672 *this->m_WriteStream << "\t % Dimension , deformed state local coordinates";
673 *this->m_WriteStream << "\n";
674
675 // print square root of Variance
676 *this->m_WriteStream << Load->m_Variance;
677 *this->m_WriteStream << "\t % Square root of the landmark variance ";
678 *this->m_WriteStream << "\n";
679 return;
680 }
681 }
682
683 void
SkipWhiteSpace()684 MetaFEMObject::SkipWhiteSpace()
685 {
686 std::string skip;
687
688 while ( this->m_ReadStream && !this->m_ReadStream->eof() && ( std::ws(*this->m_ReadStream).peek() ) == '%' )
689 {
690 std::getline(*this->m_ReadStream, skip);
691 }
692 }
693
IsClassNamePresent(std::string c_string)694 bool MetaFEMObject::IsClassNamePresent(std::string c_string)
695 {
696 ClassNameListType::const_iterator it = this->m_ClassNameList.begin();
697 while(it != this->m_ClassNameList.end())
698 {
699 if((*it) == c_string)
700 return true;
701 ++it;
702 }
703 return false;
704 }
705
M_Read_Node()706 bool MetaFEMObject::M_Read_Node()
707 {
708 unsigned int n;
709 float coor[3];
710 /**
711 * First call the parent's read function
712 */
713 int GN = this->ReadGlobalNumber();
714
715 if(GN == -1)
716 {
717 METAIO_STREAM::cout << "Error reading Global Number" << METAIO_STREAM::endl;
718 return false;
719 }
720 /*
721 * Read and set node coordinates
722 */
723 // read dimensions
724 this->SkipWhiteSpace(); *this->m_ReadStream >> n;
725 if ( !this->m_ReadStream)
726 {
727 METAIO_STREAM::cout << "Error reading Node dimensions" << METAIO_STREAM::endl;
728 return false;
729 }
730 FEMObjectNode *node = new FEMObjectNode(n);
731 node->m_GN = GN;
732
733 this->SkipWhiteSpace();
734 for (unsigned int i = 0; i< n; i++)
735 {
736 *this->m_ReadStream >> coor[i];
737 if ( !this->m_ReadStream )
738 {
739 METAIO_STREAM::cout << "Error reading Node coordinates" << METAIO_STREAM::endl;
740 return false;
741 }
742 node->m_X[i] = coor[i];
743 }
744 this->m_NodeList.push_back(node);
745
746 return true;
747 }
748
M_Read_Material(std::string material_name)749 bool MetaFEMObject::M_Read_Material(std::string material_name)
750 {
751 /**
752 * First call the parent's read function
753 */
754 int GN = this->ReadGlobalNumber();
755
756 if(GN == -1)
757 {
758 METAIO_STREAM::cout << "Error reading Global Number" << METAIO_STREAM::endl;
759 return false;
760 }
761 /*
762 * Read material properties
763 */
764 double d;
765
766 std::streampos l(0);
767 char buf[256];
768 std::string s;
769 std::string::size_type b, e;
770
771 // clear the data already inside the object
772 double E = 0.0; double A = 0.0; double I = 0.0;
773 double nu = 0.0; double h = 1.0; double RhoC = 1.0;
774
775 /*
776 * Next we read any known constant from stream. This allows a user to
777 * specify only constants which are actually required by elements in
778 * a system. This makes creating input files a bit easier.
779 */
780 while ( this->m_ReadStream )
781 {
782 l = this->m_ReadStream->tellg(); // remember the stream position
783 this->SkipWhiteSpace(); // skip comments and whitespaces
784
785 /**
786 * All Constants are in the following format:
787 * constant_name : value
788 */
789 this->m_ReadStream->getline(buf, 256, ':'); // read up to 256 characters until ':' is
790 // reached. we read and discard the ':'
791 if ( !this->m_ReadStream )
792 {
793 METAIO_STREAM::cout << "Error reading Material properties" << METAIO_STREAM::endl;
794 return false;
795 } // no : was found
796 s = std::string(buf);
797
798 // Get rid of the whitespaces in front of and the back of token
799 b = s.find_first_not_of(whitespaces); // end
800 // of
801 // whitespaces
802 // in
803 // the
804 // beginning
805 if ( ( e = s.find_first_of(whitespaces, b) ) == std::string::npos ) //
806 // beginning
807 // of
808 // whitespaces
809 // at
810 // the
811 // end
812 {
813 e = s.size();
814 }
815
816 /*
817 * s now contains just the name of the constant.
818 * The value is ready to be read next from the stream
819 */
820 s = s.substr(b, e - b);
821
822 if ( s == "E" )
823 {
824 *this->m_ReadStream >> d;
825 if ( !this->m_ReadStream )
826 {
827 METAIO_STREAM::cout << "Error reading Material E property" << METAIO_STREAM::endl;
828 return false;
829 }
830 E = d;
831 continue;
832 }
833
834 if ( s == "A" )
835 {
836 *this->m_ReadStream >> d;
837 if ( !this->m_ReadStream )
838 {
839 METAIO_STREAM::cout << "Error reading Material A property" << METAIO_STREAM::endl;
840 return false;
841 }
842 A = d;
843 continue;
844 }
845
846 if ( s == "I" )
847 {
848 this->SkipWhiteSpace();
849 *this->m_ReadStream >> d;
850 if ( !this->m_ReadStream )
851 {
852 METAIO_STREAM::cout << "Error reading Material I property" << METAIO_STREAM::endl;
853 return false;
854 }
855 I = d;
856 continue;
857 }
858
859 if ( s == "nu" )
860 {
861 this->SkipWhiteSpace();
862 *this->m_ReadStream >> d;
863 if ( !this->m_ReadStream )
864 {
865 METAIO_STREAM::cout << "Error reading Material nu property" << METAIO_STREAM::endl;
866 return false;
867 }
868 nu = d;
869 continue;
870 }
871
872 if ( s == "h" )
873 {
874 this->SkipWhiteSpace();
875 *this->m_ReadStream >> d;
876 if ( !this->m_ReadStream )
877 {
878 METAIO_STREAM::cout << "Error reading Material h property" << METAIO_STREAM::endl;
879 return false;
880 }
881 h = d;
882 continue;
883 }
884
885 if ( s == "RhoC" )
886 {
887 this->SkipWhiteSpace();
888 *this->m_ReadStream >> d;
889 if ( !this->m_ReadStream )
890 {
891 METAIO_STREAM::cout << "Error reading Material RhoC property" << METAIO_STREAM::endl;
892 return false;
893 }
894 RhoC = d;
895 continue;
896 }
897 if ( s == "END" )
898 {
899 // End of constants in material definition
900 // store all the material definitions
901 FEMObjectMaterial *material = new FEMObjectMaterial();
902 strcpy(material->m_MaterialName, material_name.c_str());
903 material->m_GN = GN;
904 material->E = E;
905 material->A = A;
906 material->I = I;
907 material->nu = nu;
908 material->h = h;
909 material->RhoC = RhoC;
910 this->m_MaterialList.push_back(material);
911 break;
912 }
913
914 /**
915 * If we got here an unknown constant was reached.
916 * We reset the stream position and set the stream error.
917 */
918 this->m_ReadStream->seekg(l);
919 this->m_ReadStream->clear(std::ios::failbit);
920 }
921
922 if ( !this->m_ReadStream )
923 {
924 METAIO_STREAM::cout << "Error reading Material properties" << METAIO_STREAM::endl;
925 return false;
926 }
927 return true;
928 }
929
M_Read_Element(std::string element_name)930 bool MetaFEMObject::M_Read_Element(std::string element_name)
931 {
932 unsigned int n, materialGN;
933 int info[2];
934 this->GetElementDimensionAndNumberOfNodes(element_name, info);
935
936 int GN = this->ReadGlobalNumber();
937
938 if(GN == -1)
939 {
940 METAIO_STREAM::cout << "Error reading Global Number" << METAIO_STREAM::endl;
941 return false;
942 }
943 /*
944 * Read and set element connectivity
945 */
946 int *NodesId = new int[info[0]];
947 for ( int p = 0; p < info[0]; p++ )
948 {
949 this->SkipWhiteSpace(); *this->m_ReadStream >> n;
950 if ( !this->m_ReadStream )
951 {
952 delete [] NodesId;
953 METAIO_STREAM::cout << "Error reading Element node numbers" << METAIO_STREAM::endl;
954 return false;
955 }
956 NodesId[p] = n;
957 }
958
959 // read material associated with the element
960 this->SkipWhiteSpace(); *this->m_ReadStream >> materialGN;
961 if ( !this->m_ReadStream )
962 {
963 delete [] NodesId;
964 METAIO_STREAM::cout << "Error reading Element global number" << METAIO_STREAM::endl;
965 return false;
966 }
967 // store the read information
968 FEMObjectElement *element = new FEMObjectElement(info[0]);
969 element->m_GN = GN;
970 for ( int p = 0; p < info[0]; p++ )
971 {
972 element->m_NodesId[p] = NodesId[p];
973 }
974 element->m_MaterialGN = materialGN;
975 element->m_NumNodes = info[0];
976 element->m_Dim = info[1];
977 strcpy(element->m_ElementName, element_name.c_str());
978
979 delete [] NodesId;
980 this->m_ElementList.push_back(element);
981 return true;
982 }
983
M_Read_Load(std::string load_name)984 bool MetaFEMObject::M_Read_Load(std::string load_name)
985 {
986 int GN;
987 int elementGN;
988 int DOF;
989 int NumRHS;
990 int NodeNumber;
991 int Dim;
992 int NumLHS;
993 int Value;
994
995 FEMObjectLoad *load = new FEMObjectLoad;
996 strcpy(load->m_LoadName, load_name.c_str());
997
998 GN = this->ReadGlobalNumber();
999 if (GN == -1)
1000 {
1001 delete load;
1002 METAIO_STREAM::cout << "Error reading Load definition - global number" << METAIO_STREAM::endl;
1003 return false;
1004 }
1005
1006 load->m_GN = GN;
1007
1008 if(load_name == "LoadBC")
1009 {
1010 /* read and set pointer to element that we're applying the load to */
1011 this->SkipWhiteSpace();
1012 *this->m_ReadStream >> elementGN;
1013 if(!this->m_ReadStream)
1014 {
1015 delete load;
1016 METAIO_STREAM::cout << "Error reading Load definition - Element Global Number" << METAIO_STREAM::endl;
1017 return false;
1018 }
1019 load->m_ElementGN = elementGN;
1020
1021 /* read the local DOF number within that element */
1022 this->SkipWhiteSpace();
1023 *this->m_ReadStream >> DOF;
1024 if(!this->m_ReadStream)
1025 {
1026 delete load;
1027 METAIO_STREAM::cout << "Error reading Load definition - Degrees of Freedom" << METAIO_STREAM::endl;
1028 return false;
1029 }
1030 load->m_DOF = DOF;
1031
1032 /* read the value to which the DOF is fixed */
1033 this->SkipWhiteSpace();
1034 *this->m_ReadStream >> NumRHS;
1035 if(!this->m_ReadStream)
1036 {
1037 delete load;
1038 METAIO_STREAM::cout << "Error reading Load definition - Number of fixed degrees of freedom" << METAIO_STREAM::endl;
1039 return false;
1040 }
1041 load->m_NumRHS = NumRHS;
1042 load->m_RHS.resize(NumRHS);
1043
1044 for (int i=0; i<NumRHS; i++)
1045 {
1046 this->SkipWhiteSpace();
1047 *this->m_ReadStream >> load->m_RHS[i];
1048 if(!this->m_ReadStream)
1049 {
1050 delete load;
1051 METAIO_STREAM::cout << "Error reading Load definition - Fixed degree of freedom" << METAIO_STREAM::endl;
1052 return false;
1053 }
1054 }
1055 }
1056 else if(load_name == "LoadNode")
1057 {
1058 /* read and set pointer to element that we're applying the load to */
1059 this->SkipWhiteSpace();
1060 *this->m_ReadStream >> elementGN;
1061 if(!this->m_ReadStream)
1062 {
1063 delete load;
1064 METAIO_STREAM::cout << "Error reading LoadNode definition - Element Global Number" << METAIO_STREAM::endl;
1065 return false;
1066 }
1067 load->m_ElementGN = elementGN;
1068
1069 /* read the value to which the DOF is fixed */
1070 this->SkipWhiteSpace();
1071 *this->m_ReadStream >> NodeNumber;
1072 if(!this->m_ReadStream)
1073 {
1074 delete load;
1075 METAIO_STREAM::cout << "Error reading LoadNode definition - Node Number" << METAIO_STREAM::endl;
1076 return false;
1077 }
1078 load->m_NodeNumber = NodeNumber;
1079
1080 this->SkipWhiteSpace();
1081 *this->m_ReadStream >> Dim;
1082 if(!this->m_ReadStream)
1083 {
1084 delete load;
1085 METAIO_STREAM::cout << "Error reading LoadNode definition - Dimension" << METAIO_STREAM::endl;
1086 return false;
1087 }
1088 load->m_Dim = Dim;
1089
1090 load->m_ForceVector.resize(Dim);
1091 for (int i=0; i<Dim; i++)
1092 {
1093 this->SkipWhiteSpace();
1094 *this->m_ReadStream >> load->m_ForceVector[i];
1095 if(!this->m_ReadStream)
1096 {
1097 delete load;
1098 METAIO_STREAM::cout << "Error reading LoadNode definition - Force Vector" << METAIO_STREAM::endl;
1099 return false;
1100 }
1101 }
1102 }
1103 else if(load_name == "LoadBCMFC")
1104 {
1105 /** read number of terms in lhs of MFC equation */
1106 this->SkipWhiteSpace();
1107 *this->m_ReadStream >> NumLHS;
1108 if(!this->m_ReadStream)
1109 {
1110 delete load;
1111 METAIO_STREAM::cout << "Error reading LoadBCMFC definition - Number of LHS terms" << METAIO_STREAM::endl;
1112 return false;
1113 }
1114
1115 load->m_NumLHS = NumLHS;
1116 for ( int i = 0; i < NumLHS; i++ )
1117 {
1118 /** read and set pointer to element that we're applying the load to */
1119 this->SkipWhiteSpace();
1120 *this->m_ReadStream >> elementGN;
1121 if(!this->m_ReadStream)
1122 {
1123 delete load;
1124 METAIO_STREAM::cout << "Error reading LoadBCMFC definition - Element Global Number" << METAIO_STREAM::endl;
1125 return false;
1126 }
1127
1128 /** read the number of dof within that element */
1129 this->SkipWhiteSpace();
1130 *this->m_ReadStream >> DOF;
1131 if(!this->m_ReadStream)
1132 {
1133 delete load;
1134 METAIO_STREAM::cout << "Error reading LoadBCMFC definition - Element Degree of Freedom" << METAIO_STREAM::endl;
1135 return false;
1136 }
1137
1138 /** read weight */
1139 this->SkipWhiteSpace();
1140 *this->m_ReadStream >> Value;
1141 if(!this->m_ReadStream)
1142 {
1143 delete load;
1144 METAIO_STREAM::cout << "Error reading LoadBCMFC definition - Weight" << METAIO_STREAM::endl;
1145 return false;
1146 }
1147
1148 /** add a new MFCTerm to the lhs */
1149 FEMObjectMFCTerm *mfcTerm = new FEMObjectMFCTerm(elementGN, DOF, static_cast<float>(Value));
1150 load->m_LHS.push_back(mfcTerm);
1151 }
1152
1153 /** read the rhs */
1154 this->SkipWhiteSpace();
1155 *this->m_ReadStream >> NumRHS;
1156 if (!this->m_ReadStream)
1157 {
1158 delete load;
1159 METAIO_STREAM::cout << "Error reading LoadBCMFC definition - Number of RHS terms" << METAIO_STREAM::endl;
1160 return false;
1161 }
1162
1163 load->m_NumRHS = NumRHS;
1164 load->m_RHS.resize(NumRHS);
1165 for (int i=0; i<NumRHS; i++)
1166 {
1167 this->SkipWhiteSpace();
1168 *this->m_ReadStream >> load->m_RHS[i];
1169 if(!this->m_ReadStream)
1170 {
1171 delete load;
1172 METAIO_STREAM::cout << "Error reading LoadBCMFC definition - RHS Term" << METAIO_STREAM::endl;
1173 return false;
1174 }
1175 }
1176 }
1177 else if(load_name == "LoadEdge")
1178 {
1179 int edgeNum, numRows, numCols;
1180
1181 /* read the global number of the element on which the load acts */
1182 this->SkipWhiteSpace();
1183 *this->m_ReadStream >> elementGN;
1184 if (!this->m_ReadStream)
1185 {
1186 delete load;
1187 METAIO_STREAM::cout << "Error reading LoadEdge definition - Element Global Number" << METAIO_STREAM::endl;
1188 return false;
1189 }
1190 load->m_ElementGN = elementGN;
1191
1192 /** ... edge number */
1193 this->SkipWhiteSpace();
1194 *this->m_ReadStream >> edgeNum;
1195 if (!this->m_ReadStream)
1196 {
1197 delete load;
1198 METAIO_STREAM::cout << "Error reading LoadEdge definition - Edge Number" << METAIO_STREAM::endl;
1199 return false;
1200 }
1201 load->m_EdgeNumber = edgeNum;
1202
1203 /** ... # of rows */
1204 this->SkipWhiteSpace();
1205 *this->m_ReadStream >> numRows;
1206 if (!this->m_ReadStream)
1207 {
1208 delete load;
1209 METAIO_STREAM::cout << "Error reading LoadEdge definition - Number of Rows" << METAIO_STREAM::endl;
1210 return false;
1211 }
1212
1213 /** ... # of cols */
1214 this->SkipWhiteSpace();
1215 *this->m_ReadStream >> numCols;
1216 if (!this->m_ReadStream)
1217 {
1218 delete load;
1219 METAIO_STREAM::cout << "Error reading LoadEdge definition - Number of Columns" << METAIO_STREAM::endl;
1220 return false;
1221 }
1222
1223 for ( int i = 0; i < numRows; i++ )
1224 {
1225 this->SkipWhiteSpace();
1226 METAIO_STL::vector<float> F(numCols);
1227 for ( int j = 0; j < numCols; j++ )
1228 {
1229 *this->m_ReadStream >> F[j];
1230 if (!this->m_ReadStream)
1231 {
1232 delete load;
1233 METAIO_STREAM::cout << "Error reading LoadEdge definition - Force Matrix" << METAIO_STREAM::endl;
1234 return false;
1235 }
1236 }
1237 this->SkipWhiteSpace();
1238 load->m_ForceMatrix.push_back(F);
1239 }
1240 }
1241 else if(load_name == "LoadGravConst")
1242 {
1243 // read in the list of elements on which the load acts
1244 this->SkipWhiteSpace();
1245 *this->m_ReadStream >> load->m_NumElements ;
1246 if (!this->m_ReadStream)
1247 {
1248 delete load;
1249 METAIO_STREAM::cout << "Error reading LoadGravConst definition - Number of Elements" << METAIO_STREAM::endl;
1250 return false;
1251 }
1252
1253 for (int i=0; i<load->m_NumElements; i++)
1254 {
1255 this->SkipWhiteSpace();
1256 *this->m_ReadStream >> elementGN ;
1257 if (!this->m_ReadStream)
1258 {
1259 delete load;
1260 METAIO_STREAM::cout << "Error reading LoadGravConst definition - Element Global Number" << METAIO_STREAM::endl;
1261 return false;
1262 }
1263 load->m_Elements.push_back(elementGN);
1264 }
1265
1266 /** first read and set the size of the vector */
1267 this->SkipWhiteSpace();
1268 *this->m_ReadStream >> load->m_Dim;
1269 if (!this->m_ReadStream)
1270 {
1271 delete load;
1272 METAIO_STREAM::cout << "Error reading LoadGravConst definition - Dimension" << METAIO_STREAM::endl;
1273 return false;
1274 }
1275
1276 float loadcomp;
1277 /** then the actual values */
1278 for (int i=0; i<load->m_Dim; i++)
1279 {
1280 this->SkipWhiteSpace();
1281 *this->m_ReadStream >> loadcomp;
1282 if (!this->m_ReadStream)
1283 {
1284 delete load;
1285 METAIO_STREAM::cout << "Error reading LoadGravConst definition - Force Vector" << METAIO_STREAM::endl;
1286 return false;
1287 }
1288 load->m_ForceVector.push_back(loadcomp);
1289 }
1290 }
1291 else if(load_name == "LoadLandmark")
1292 {
1293 this->SkipWhiteSpace();
1294 int n1, n2;
1295
1296 // read the dimensions of the undeformed point and set the size of the point
1297 // accordingly
1298 this->SkipWhiteSpace();
1299 *this->m_ReadStream >> n1; if ( !this->m_ReadStream ) { return false; }
1300 load->m_Undeformed.resize(n1);
1301 for (int i=0; i<n1; i++)
1302 {
1303 this->SkipWhiteSpace();
1304 *this->m_ReadStream >> load->m_Undeformed[i];
1305 std::cout << " " << load->m_Undeformed[i] << std::endl;
1306 if(!this->m_ReadStream)
1307 {
1308 delete load;
1309 METAIO_STREAM::cout << "Error reading Loadlandmark definition - Undeformed point" << METAIO_STREAM::endl;
1310 return false;
1311 }
1312 }
1313
1314 // Read the dimensions of the deformed point and set the size of the point
1315 // accordingly
1316 this->SkipWhiteSpace();
1317 *this->m_ReadStream >> n2; if ( !this->m_ReadStream ) { return false; }
1318 load->m_Deformed.resize(n2);
1319 for (int i=0; i<n2; i++)
1320 {
1321 this->SkipWhiteSpace();
1322 *this->m_ReadStream >> load->m_Deformed[i];
1323 std::cout << " " << load->m_Deformed[i] << std::endl;
1324 if(!this->m_ReadStream)
1325 {
1326 delete load;
1327 METAIO_STREAM::cout << "Error reading Loadlandmark definition - Undeformed point" << METAIO_STREAM::endl;
1328 return false;
1329 }
1330 }
1331
1332 // Verify that the undeformed and deformed points are of the same size.
1333 if ( n1 != n2 )
1334 {
1335 delete load;
1336 METAIO_STREAM::cout << "Error reading Loadlandmark definition - Undeformed point and deformed point should have same dimension" << METAIO_STREAM::endl;
1337 return false;
1338 }
1339 // read the square root of the m_Variance associated with this landmark
1340 this->SkipWhiteSpace();
1341 *this->m_ReadStream >> load->m_Variance;
1342 }
1343 if(!this->m_ReadStream)
1344 {
1345 delete load;
1346 METAIO_STREAM::cout << "Error reading Load definition" << METAIO_STREAM::endl;
1347 return false;
1348 }
1349 this->m_LoadList.push_back(load);
1350 return true;
1351 }
1352
GetElementDimensionAndNumberOfNodes(std::string c_string,int info[2])1353 int* MetaFEMObject::GetElementDimensionAndNumberOfNodes(std::string c_string, int info[2])
1354 {
1355 if((c_string == "Element2DC0LinearLineStress") || (c_string == "Element2DC1Beam"))
1356 {
1357 info[0] = 2;
1358 info[1] = 2;
1359 }
1360
1361 if((c_string == "Element2DC0LinearTriangularMembrane") ||
1362 (c_string == "Element2DC0LinearTriangularStrain") ||
1363 (c_string == "Element2DC0LinearTriangularStress"))
1364 {
1365 info[0] = 3;
1366 info[1] = 2;
1367 }
1368
1369 if((c_string == "Element2DC0LinearQuadrilateralMembrane") ||
1370 (c_string == "Element2DC0LinearQuadrilateralStrain") ||
1371 (c_string == "Element2DC0LinearQuadrilateralStress"))
1372 {
1373 info[0] = 4;
1374 info[1] = 2;
1375 }
1376
1377
1378 if((c_string == "Element2DC0QuadraticTriangularStrain") ||
1379 (c_string == "Element2DC0QuadraticTriangularStress"))
1380 {
1381 info[0] = 6;
1382 info[1] = 2;
1383 }
1384
1385 if((c_string == "Element3DC0LinearHexahedronMembrane") ||
1386 (c_string == "Element3DC0LinearHexahedronStrain"))
1387 {
1388 info[0] = 8;
1389 info[1] = 3;
1390 }
1391
1392 if((c_string == "Element3DC0LinearTetrahedronMembrane") ||
1393 (c_string == "Element3DC0LinearTetrahedronStrain"))
1394 {
1395 info[0] = 4;
1396 info[1] = 3;
1397 }
1398
1399 return info;
1400 }
ReadGlobalNumber()1401 int MetaFEMObject::ReadGlobalNumber()
1402 {
1403 int n;
1404
1405 /** Read and set the global object number */
1406 this->SkipWhiteSpace();
1407 *this->m_ReadStream >> n;
1408 if(this->m_ReadStream)
1409 return n;
1410 else
1411 return -1;
1412 }
1413
1414 // string containing all whitespace characters
1415 const std::string
1416 MetaFEMObject
1417 :: whitespaces = " \t\n\r";
1418
1419 #if (METAIO_USE_NAMESPACE)
1420 }
1421 #endif
1422