/*********************************************************************** * * ELMER, A Computational Fluid Dynamics Program. * * Copyright 1st April 1995 - , Center for Scientific Computing, * Finland. * * All rights reserved. No part of this program may be used, * reproduced or transmitted in any form or by any means * without the written permission of CSC. * * Address: Center for Scientific Computing * Tietotie 6, P.O. BOX 405 * 02101 Espoo, Finland * Tel. +358 0 457 2001 * Telefax: +358 0 457 2302 * EMail: Jari.Jarvinen@csc.fi ************************************************************************/ /*********************************************************************** Program: ELMER Front Module: ecif_inputFidap.cpp Language: C++ Date: 27.04.99 Version: 1.00 Author(s): Martti Verho Revisions: Abstract: Implementation ************************************************************************/ #include "ecif_body.h" #include "ecif_body2D.h" #include "ecif_body3D.h" #include "ecif_bodyElement.h" #include "ecif_bodyElement2D.h" #include "ecif_bodyElement3D.h" #include "ecif_control.h" #include "ecif_geometry.h" #include "ecif_inputFidap.h" #include "ecif_mesh.h" #include "ecif_model.h" #include "ecif_timer.h" #include "ecif_userinterface.h" extern char read_buffer[]; // Indices to the special reorder table for element types const int fidapReorderTableIndex_203 = 0; const int fidapReorderTableIndex_510 = 1; const int fidapReorderTableIndex_718 = 2; const int fidapReorderTableIndex_808 = 3; const int fidapReorderTableIndex_809 = 4; const int fidapReorderTableIndex_827 = 5; // We have to reorder Fidap nodes to match Elmer numbering const int fidapReorderTable[][MAX_NOF_NODES] = { {0,2,1}, // Line 203 {0,2,5,9,1,4,3,6,7,8}, // Tetra 510 // NOTE: 718 is delivered to Elmer as 706, but here is the complete reorder! (ok?) {0,2,5,12,14,17,1,4,3,6,7,8,10,11,9,13,16,15}, // Wedge 718 {0,1,3,2,4,5,7,6}, // Brick 808 {0,1,3,2,4,5,7,6, 8}, // Brick 809 (1 in center) // NOTE: 827 is delivered to Elmer as 820, but here is the complete reorder! {0,2,8,6,18,20,26,24,1,5,7,3,9,11,17,15,19,23,25,21,10,14,16,12,4,22,13} // Brick 827 }; InputFidap::InputFidap(enum ecif_modelDimension m_dim, ifstream& in_file, char* in_filename): Input(m_dim, in_file, in_filename) { for (short i = 0; i < 1 + MAX_NOF_ELEM_CODES; i++) elementCodeCounters[i] = 0; nofElementGroups = 0; dataDimension = ECIF_ND; } InputFidap::~InputFidap() { } Body* InputFidap::findExistingBody(FidapElementGroupInfo* gi) { int index = 0; while (true) { Body* body = model->getBody(index++); if (body==NULL) break; if ( LibFront::ncEqual((char*)body->getName(), gi->entityName) ) { return body; } } return NULL; } BodyElement* InputFidap::findExistingBoundary(FidapElementGroupInfo* gi) { int index = 0; while (true) { BodyElement* be = model->getBoundary(index++); if (be==NULL) break; if ( LibFront::ncEqual((char*)be->getName(), gi->entityName) ) { return be; } } return NULL; } Rc InputFidap::findMaxExternalIds(int& max_ext_nd_id, int& max_ext_el_id) { enum fidapKeywordId keyword = FIDAP_NO_ID; bool count_only = true; while ( !infile.eof() ) { // Read new line readFileLine(infile, read_buffer); // Check if it was an interesting keyword line keyword = next_is_keyword_line(read_buffer); if ( keyword == FIDAP_NO_ID) continue; // Read element group // ================== if ( keyword == FIDAP_ELEMENT_GROUPS_ID ) { int element_group_counter = 0; while ( element_group_counter < nofElementGroups && !infile.eof() ) { element_group_counter++; // Read one group FidapElementGroupInfo gi; // Check that group is ok if ( !readElementGroupInfo(gi) ) return ECIF_ERROR; if ( !setElementType(gi) ) return ECIF_ERROR; // Read group elements // =================== int element_counter = 0; readElements(gi, element_counter, max_ext_el_id, count_only); } } // Read nodes // ========== else if ( keyword == FIDAP_NODAL_COORDINATES_ID ) { int node_counter = 0; readNodes(nofNodes, node_counter, max_ext_nd_id, count_only); } } // while !eof return ECIF_OK; } // NOTE: Fidap mesh dimension is read from the header and // stored in the class attribute 'dimension' // enum ecif_modelDimension InputFidap::findMeshModelDimension() { ecif_modelDimension box_dim = ECIF_ND; ecif_modelDimension dim = ECIF_ND; #if 0 RangeVector rv; bool x_is_const = false; bool y_is_const = false; bool z_is_const = false; meshBoundingBox.getRangeVector(rv); // Find model dimension by boundary box if ( isEqual(rv[0], rv[1]) ) { x_is_const = true; } if ( isEqual(rv[2], rv[3]) ) { y_is_const = true; } if ( isEqual(rv[4], rv[5]) ) { z_is_const = true; } if ( !(x_is_const || y_is_const || z_is_const) ) { box_dim = ECIF_3D; } else if ( !( x_is_const || y_is_const) && z_is_const ) { box_dim = ECIF_2D; } else if ( !(x_is_const && y_is_const && z_is_const) ) { box_dim = ECIF_1D; } else { box_dim = ECIF_ND; } #endif box_dim = modelDimension; // Find actual dimension dim = findModelDimension(box_dim); return dim; } enum fidapKeywordId InputFidap::next_is_keyword_line(char* line_buffer) { if ( 0 == strncmp("NODAL COORDINATES", line_buffer, 17) ) return FIDAP_NODAL_COORDINATES_ID; else if ( 0 == strncmp("ELEMENT GROUPS", line_buffer, 14) ) return FIDAP_ELEMENT_GROUPS_ID; else if ( 0 == strncmp("VERSION", line_buffer, 7) ) return FIDAP_VERSION_ID; else return FIDAP_NO_ID; } Rc InputFidap::readElementGroup(int& body_counter, int& bndr_counter, int& element_counter, int& max_ext_elem_id, bool count_only) { int element_group_counter = 0; while ( element_group_counter < nofElementGroups && !infile.eof() ) { element_group_counter++; // Read one group FidapElementGroupInfo gi; // Check that group is ok if ( !readElementGroupInfo(gi) ) return ECIF_ERROR; if ( !setElementType(gi) ) return ECIF_ERROR; gi.parentTag = NO_INDEX; // Check elment type bulk/boundary model->checkMeshElementType(gi.elementType, gi.isBulk, gi.isBoundary, gi.isEdge); // Bulk elements --> add elements to a body // ============= if (gi.isBulk) { nofBulkElements += gi.nofElements; // Try to find an existing body // ---------------------------- Body* body = findExistingBody(&gi); // Create a new body // ----------------- if (body == NULL) { int int_id = ++body_counter; // Body internal id! int ext_id = gi.id; // Body external id! if ( modelDimension == ECIF_2D ) { body = new Body2D(MESH_BODY, int_id, ext_id, 0, NULL); } else if ( modelDimension == ECIF_3D ) { body = new Body3D(MESH_BODY, int_id, ext_id, 0, NULL); } if (body != NULL) { model->addBody(body); body->setName(gi.entityName); } } // Ok, home found/created for the bulk elements! if ( body != NULL ) { gi.parentTag = body->Tag(); } // Boundary elements --> add elements to a boundary // ================= } else if (gi.isBoundary) { nofBoundaryElements += gi.nofElements; // Try to find an existing boundary // -------------------------------- BodyElement* be = findExistingBoundary(&gi); // Create a new boundary // --------------------- if ( be == NULL ) { int int_id = ++bndr_counter; // Internal id! // Create a proper boundary (body element) if ( modelDimension == ECIF_2D ) { be = new BodyElement2D(int_id, NO_INDEX, NO_INDEX, gi.nofElements); } else if ( modelDimension == ECIF_3D ) { be = new BodyElement3D(int_id, NO_INDEX, NO_INDEX, gi.nofElements); } if (be != NULL) { model->addBodyElement(be); be->setName(gi.entityName); } } // Ok, home found/created for the boundary elements! if ( be != NULL ) { gi.parentTag = be->Tag(); } } // Read group elements // =================== if ( ECIF_OK != readElements(gi, element_counter, max_ext_elem_id, count_only) ) { return ECIF_ERROR; } } return ECIF_OK; } // Read element group info // Returns true if read was succesful // bool InputFidap::readElementGroupInfo(FidapElementGroupInfo& gi) { static char key_buffer[80]; static char value_buffer[80]; eat_white(infile); // First group info line readFileLine(infile, read_buffer); strstream strm1; strm1 << read_buffer; while( !strm1.eof() ) { ws(strm1); strm1.getline(key_buffer, 80, ':'); if ( 0 == strncmp(key_buffer, "GROUP", 5) ) strm1 >> gi.id; if ( 0 == strncmp(key_buffer, "ELEMENTS", 8) ) strm1 >> gi.nofElements; if ( 0 == strncmp(key_buffer, "NODES", 5) ) strm1 >> gi.nofNodes; if ( 0 == strncmp(key_buffer, "GEOMETRY", 8) ) strm1 >> gi.fidap_geometry; if ( 0 == strncmp(key_buffer, "TYPE", 4) ) strm1 >> gi.fidap_type; } // Second group info line readFileLine(infile, read_buffer); strstream strm2; strm2 << read_buffer; while( !strm2.eof() ) { ws(strm2); strm2.getline(key_buffer, 80, ':'); if ( 0 == strncmp(key_buffer, "ENTITY NAME", 11) ) strm2.getline(gi.entityName, 80 , '\n'); } LibFront::trim(gi.entityName); return true; } Rc InputFidap::readElements(FidapElementGroupInfo& gi, int& element_counter, int& max_ext_id, bool count_only) { Rc rc; int ext_el_id, ext_nd_id; static int ext_node_ids[MAX_NOF_NODES]; static int ext_node_ids2[MAX_NOF_NODES]; // For reordering! meshElementCode elem_code = model->convertElementType(gi.elementType); if ( elem_code == MEC_000 ) { return ECIF_ERROR; } short nof_nodes = MeshElementDesc[elem_code][DESC_NOF_NODES]; short i, j; eat_white(infile); int group_elem_counter = 0; while ( group_elem_counter < gi.nofElements && !infile.eof() ) { if ( !count_only ) { theControlCenter->update(element_counter, GUI_UPDATE_INTERVAL); } int node_count = 0; while ( node_count < nof_nodes ) { // Element id if ( node_count == 0 ) { infile >> ext_el_id; } infile >> ext_nd_id; ext_node_ids[node_count] = ext_node_ids2[node_count] = ext_nd_id; node_count++; } if ( max_ext_id < ext_el_id ) { max_ext_id = ext_el_id; } group_elem_counter++; element_counter++; if ( count_only ) continue; // Apply possible preset reordering vector if ( gi.fidap_reorder != NULL ) { for (i = 0; i < nof_nodes; i++) { ext_node_ids[i] = ext_node_ids2[ gi.fidap_reorder[i] ]; } // Or apply possible generic reorder for non-ok element types } else { switch (gi.elementType) { // These are ok case 303: case 304: case 404: case 405: case 504: case 505: case 706: case 707: break; // Others have a generic reordering default: for (i = 0, j = 0; j < nof_nodes - 1; i++, j += 2) { ext_node_ids[i] = ext_node_ids2[j]; } for (j = 1; j < nof_nodes; i++, j += 2) { ext_node_ids[i] = ext_node_ids2[j]; } break; } } elementCodeCounters[elem_code]++; // Add element to the model rc = model->addMeshInputElement(gi.elementType, ext_el_id, gi.parentTag, gi.id, ext_node_ids); if (gi.isBulk) { nofInputBulkElements++; } else if (gi.isBoundary) { nofInputBoundaryElements++; } else if (gi.isEdge) { nofInputEdgeElements++; } else { nofInputVertexElements++; } if (rc != ECIF_OK) { return rc; } } // while not eof() return ECIF_OK; } bool InputFidap::readNumericKeywordValue(char* line_buffer, char* keyword, int& value) { const int buffer_len = 80; static char buffer[buffer_len]; if ( !readKeywordValue(line_buffer, keyword, buffer, buffer_len) ) return false; value = atol(buffer); return true; } bool InputFidap::readStringKeywordValue(char* line_buffer, char* keyword, char* value_buffer) { const int buffer_len = 80; static char buffer[buffer_len]; if ( !readKeywordValue(line_buffer, keyword, buffer, buffer_len) ) return false; strcpy(buffer, value_buffer); return true; } bool InputFidap::readKeyword(char* line_buffer, char* buffer, int buffer_len) { strstream strm; strm << line_buffer << ends; if ( strm.eof() ) return false; strm.getline(buffer, buffer_len, ':'); // Trim trailing blanks int len = strlen(buffer); int pos = len - 1; while ( pos >= 0 && buffer[pos] == ' ' ) pos--; buffer[1 + pos] = '\0'; return true; } bool InputFidap::readKeywordValue(char* line_buffer, char* keyword, char* buffer, int buffer_len) { int strm_buf_len = strlen(line_buffer); char* strm_buf = new char[1 + strm_buf_len]; strstream strm(strm_buf, 1 + strm_buf_len, ios::app); strm << line_buffer << ends; // Work buffer char* work_buf = new char[1 + strm_buf_len]; char* tmp = work_buf; bool found = false; // Try to read keyword value while ( !strm.eof() ) { // Read next keyword strm.getline(tmp, strm_buf_len, ':'); // Skip leading blanks LibFront::trimLeft(tmp); // If this was searched keyword if ( LibFront::ncEqualPartial(tmp, keyword) ) { found = true; // Jump over the searched keyword tmp += strlen(keyword); // Make a new stream from what is left reset(strm); strm << tmp << ends; break; } } if (!found) return false; // Pick until next keyword to get the keyword value // for the previous keyword! strm.getline(tmp, strm_buf_len, ':'); // Skip backword the next keyword (until first blank) int len = strlen(tmp); while (len > 0 && tmp[len - 1] != ' ' ) len--; tmp[len] = '\0'; // Trim blanks LibFront::trim(tmp); // Check that value fits to the final buffer if ( strlen(tmp) <= buffer_len ) { strcpy(buffer, tmp); found = true; } else { found = false; } delete[] strm_buf; delete[] work_buf; return found; } // NOTE: Model dimension is stored in this stage in the class // attribute 'dimension' which is read from the mesh file header // Rc InputFidap::readMeshData(int& element_counter, int& node_counter) { enum fidapKeywordId keyword = FIDAP_NO_ID; bool count_only = false; UserInterface* gui = theControlCenter->getGui(); int body_counter = 0; int bndr_counter = 0; int max_ext_el_id = -1; int max_ext_nd_id = -1; while ( !infile.eof() ) { // Read new line readFileLine(infile, read_buffer); // Check if it was an interesting keyword line keyword = next_is_keyword_line(read_buffer); if ( keyword == FIDAP_NO_ID) continue; // Read elements (group) // ===================== if ( keyword == FIDAP_ELEMENT_GROUPS_ID ) { if ( ECIF_OK != readElementGroup(body_counter, bndr_counter, element_counter, max_ext_el_id, count_only) ) { return ECIF_ERROR; } } // Read nodes // ========== else if ( keyword == FIDAP_NODAL_COORDINATES_ID ) { if ( ECIF_OK != readNodes(nofNodes, node_counter, max_ext_nd_id, count_only) ) { return ECIF_ERROR; } model->setMeshNodes(); } } // while !eof return ECIF_OK; } // Create mesh geometry from Fidap neutral input file // // --Data should be in the following order in the input file: // Nodes, Elements, Element groups // // --Those bulk and boundary elements which are not listed in groups, are processed later // and the bodies/boundaries based on them are then created // // --Model dimension is known when nodes are read // bool InputFidap::readMeshGeometry() { int dimension; //---Header: Read dimension, nof nodes and elements // if (!readMeshHeader(dimension)) { modelDimension = ECIF_ND; infile.close(); return false; } if ( nofElements == 0 || nofNodes == 0 ) { modelDimension = ECIF_ND; return false; } if ( dimension == 1 ) { dataDimension = ECIF_1D; } else if ( dimension == 2 ) { dataDimension = ECIF_2D; } else if (dimension == 3 ) { dataDimension = ECIF_3D; } else { dataDimension = ECIF_ND; return false; } // Check input dimension, it cannot be smaller than // data dimension! // 1D not actually supported! if ( inputDimension == ECIF_1D && dataDimension == ECIF_2D ) { inputDimension = ECIF_2D; } if ( inputDimension == ECIF_2D && dataDimension == ECIF_3D ) { inputDimension = ECIF_3D; } modelDimension = inputDimension; model->setModelDimension(modelDimension); findMaxExternalIds(maxExternalNodeId, maxExternalElementId); //---Allocate mesh tables // NOTE: This is total of bulk+bndr elements given explicitely in the file model->allocateMeshInputElements(nofElements, maxExternalElementId); model->allocateMeshNodes(nofNodes, maxExternalNodeId); //---Read nodes and element groups Rc rc; infile.clear(); infile.seekg(0); int element_counter = 0; int node_counter = 0; rc = readMeshData(element_counter, node_counter); if (rc != ECIF_OK) { modelDimension = ECIF_ND; return false; } model->setModelDimension(modelDimension); if (modelDimension == ECIF_ND) { return false; } return true; } // All Fidap neutral files have similar header section // // NOTE: In header Nof-elements is sum of listed bulk and bndr elements! // bool InputFidap::readMeshHeader(int& dimension) { readFileLine(infile, read_buffer); if ( 0 != strncmp("** FIDAP NEUTRAL FILE", read_buffer, 21) ) { return false; } readFileLine(infile, read_buffer, 5); strstream strm; strm << read_buffer; if ( !(strm >> nofNodes >> nofElements >> nofElementGroups >> dimension) ) { return false; } return true; } Rc InputFidap::readNodes(int nof_nodes, int& node_counter, int& max_ext_id, bool count_only) { Rc rc; int ext_nd_id; static Point3 point; // Init point buffer point[0] = point[1] = point[2] = 0.0; bool is_3D = false; if ( dataDimension == ECIF_3D ) { is_3D = true; } // Read nof-nodes entries while ( !infile.eof() && node_counter < nof_nodes ) { theControlCenter->update(node_counter, GUI_UPDATE_INTERVAL); infile >> ext_nd_id; infile >> point[0] >> point[1]; if (is_3D) { infile >> point[2]; } node_counter++; if ( max_ext_id < ext_nd_id ) { max_ext_id = ext_nd_id; } // Update only counter etc. // if ( count_only ) { meshBoundingBox.extendByPoint(point); // Add node to the model // } else { rc = model->addMeshNode(node_counter - 1, ext_nd_id, point); } } // while not eof() // Model dimension can now be set! // if ( !count_only ) { modelDimension = findMeshModelDimension(); } return ECIF_OK; } // Find ELMER element code for an Fidad element bool InputFidap::setElementType(FidapElementGroupInfo& gi) { int base = 200; gi.fidap_reorder = NULL; switch (gi.fidap_geometry) { case 0: if ( dataDimension == ECIF_3D ) base = 400; break; case 1: base = 400; break; case 2: base = 300; break; case 3: base = 800; break; case 4: base = 700; break; case 5: base = 500; break; default: return false; } int type = base + gi.nofNodes; // Types needing specific reordering or mapping to // "lower" level types if (type == 203) { gi.fidap_reorder = fidapReorderTable[fidapReorderTableIndex_203]; } else if (type == 510) { gi.fidap_reorder = fidapReorderTable[fidapReorderTableIndex_510]; } else if (type >= 718 && type < 800) { if (type == 718) gi.fidap_reorder = fidapReorderTable[fidapReorderTableIndex_718]; type = 706; } else if (type == 808) { gi.fidap_reorder = fidapReorderTable[fidapReorderTableIndex_808]; } else if (type == 809) { gi.fidap_reorder = fidapReorderTable[fidapReorderTableIndex_809]; } else if (type >= 820) { if (type == 827) gi.fidap_reorder = fidapReorderTable[fidapReorderTableIndex_827]; type = 820; } gi.elementType = type; return true; }