diff --git a/NonLinearSolver/modelReduction/DeepMaterialNetworks.cpp b/NonLinearSolver/modelReduction/DeepMaterialNetworks.cpp index 07fbc4265ce7d72423706d7cd8549eac24c5aecd..27d6a83bf22fe7872d8c83f79f2b53fb5cbc2de4 100644 --- a/NonLinearSolver/modelReduction/DeepMaterialNetworks.cpp +++ b/NonLinearSolver/modelReduction/DeepMaterialNetworks.cpp @@ -722,7 +722,8 @@ void BimaterialHomogenization::computeBimaterial(const fullMatrix<double>& C1, c fullMatrix<double>* DCeffDnormal1,fullMatrix<double>* DCeffDnormal2, fullMatrix<double>* DCeffDnormal3) const { double tol = 1e-10; - if (C1.norm() < tol*C2.norm()) + double tolF = 1e-6; + if ((C1.norm() < tol*C2.norm()) || (f1 < tolF*f2)) { // C1 is void Ceff = C2; @@ -743,7 +744,7 @@ void BimaterialHomogenization::computeBimaterial(const fullMatrix<double>& C1, c fullMatrixOperation::allocateAndMakeZero(*DCeffDnormal3,Ceff.size1(),Ceff.size2()); }; } - else if (C2.norm() < tol*C1.norm()) + else if ((C2.norm() < tol*C1.norm()) || (f2 < tolF*f1)) { Ceff = C1; Ceff.scale(f1); diff --git a/NonLinearSolver/modelReduction/Tree.cpp b/NonLinearSolver/modelReduction/Tree.cpp index ed8d0a0e8ed4f125dac7a22c8804d6a51b0f1c53..208855b6941784e020ac16b71e83d85fde1988cc 100644 --- a/NonLinearSolver/modelReduction/Tree.cpp +++ b/NonLinearSolver/modelReduction/Tree.cpp @@ -33,7 +33,7 @@ TreeNode::~TreeNode() void TreeNode::printData() const { - printf("Node: depth %d- loc %d -mat index %d- child order %d - weight = %e - nb of direction vars %d",depth,location,matIndex,childOrder,weight,direction.size()); + printf("Node: depth %d- loc %d -mat index %d- child order %d - weight = %e - nb of direction vars %ld",depth,location,matIndex,childOrder,weight,direction.size()); for (int i=0; i< direction.size(); i++) { printf(" %e",direction[i]); @@ -49,7 +49,7 @@ void TreeNode::saveToFile(FILE* f) const { parentloc = parent->location; } - fprintf(f,"%d %d %d %d %d %d %d",depth,location, matIndex, childOrder,parentloc,childs.size(), direction.size()); + fprintf(f,"%d %d %d %d %d %ld %ld",depth,location, matIndex, childOrder,parentloc,childs.size(), direction.size()); // if (childs.size() > 0) { @@ -163,7 +163,6 @@ void Tree::createSpecialBinaryTree(int numNodes, int numDirVars, const std::stri { clear(); _root = new TreeNode(0,0,NULL,0); - _nodeMap[std::pair<int,int>(0,0)]=_root; TreeNode* np = _root; int depth = 1; @@ -176,7 +175,6 @@ void Tree::createSpecialBinaryTree(int numNodes, int numDirVars, const std::stri np->childs.resize(2); np->direction.resize(numDirVars); np->childs[0] = new TreeNode(depth,0,np,0,affunc); - _nodeMap[std::pair<int,int>(depth,0)]=np->childs[0]; np->childs[0]->matIndex = matStart; if (numAllocatedNodes >= numNodes) @@ -195,7 +193,6 @@ void Tree::createSpecialBinaryTree(int numNodes, int numDirVars, const std::stri { np->childs[1] = new TreeNode(depth,1,np,1); } - _nodeMap[std::pair<int,int>(depth,1)]=np->childs[1]; if (matStart == 0) matStart = 1; else matStart= 0; depth++; @@ -206,9 +203,7 @@ void Tree::createSpecialBinaryTree(int numNodes, int numDirVars, const std::stri void Tree::createSpecialTripleTree(int numNodes, int numDirVars, const std::string affunc) { clear(); - _root = new TreeNode(0,0,NULL,0); - _nodeMap[std::pair<int,int>(0,0)]=_root; - + _root = new TreeNode(0,0,NULL,0); TreeNode* np = _root; int depth = 1; int numAllocatedNodes = 1; @@ -227,16 +222,13 @@ void Tree::createSpecialTripleTree(int numNodes, int numDirVars, const std::stri np->direction.resize(numDirVars); } np->childs[0] = new TreeNode(depth,0,np,0,affunc); - _nodeMap[std::pair<int,int>(depth,0)]=np->childs[0]; np->childs[0]->matIndex = 0; np->childs[1] = new TreeNode(depth,1,np,1,affunc); - _nodeMap[std::pair<int,int>(depth,1)]=np->childs[1]; np->childs[1]->matIndex = 1; if (numAllocatedNodes < numNodes) { np->childs[2] = new TreeNode(depth,2,np,2); - _nodeMap[std::pair<int,int>(depth,2)]=np->childs[2]; np = np->childs[2]; } @@ -328,7 +320,7 @@ void Tree::createRandomTreeSameDepthForNodes(int nbLeaves, int numChildMax, int auto printVector = [](const std::vector<int>& vv) { - printf("%d:",vv.size()); + printf("%ld:",vv.size()); for (int i=0; i< vv.size(); i++) { printf("%d ",vv[i]); @@ -370,7 +362,6 @@ void Tree::createRandomTreeSameDepthForNodes(int nbLeaves, int numChildMax, int std::vector<TreeNode*> listPrev, listCur; _root = new TreeNode(0,0,NULL,0); - _nodeMap[std::pair<int,int>(0,0)]=_root; listPrev.push_back(_root); int depth = 1; @@ -400,7 +391,6 @@ void Tree::createRandomTreeSameDepthForNodes(int nbLeaves, int numChildMax, int np->childs[j] = new TreeNode(depth,inode,np,j); listCur.push_back(np->childs[j]); } - _nodeMap[std::pair<int,int>(depth,inode)]=np->childs[j]; inode++; } } @@ -542,10 +532,7 @@ void Tree::createRandomTree(int nbLeaves, int numChildMax, int leafPerNode, int std::vector<TreeNode*> listPrev, listCur; std::vector<int> listCurLeaves, listPrevLeaves; - _root = new TreeNode(0,0,NULL,0); - _nodeMap[std::pair<int,int>(0,0)]=_root; - listPrev.push_back(_root); listPrevLeaves.push_back(nbLeaves); @@ -598,7 +585,6 @@ void Tree::createRandomTree(int nbLeaves, int numChildMax, int leafPerNode, int listCurLeaves.push_back(partIndexes[j]); listCur.push_back(np->childs[j]); } - _nodeMap[std::pair<int,int>(depth,inode)]=np->childs[j]; inode++; } } @@ -667,7 +653,6 @@ void Tree::createMixedTree(int maxDepth, int numChildMax, int leafPerNode, int n if (depth == 0) { _root = new TreeNode(0,0,NULL,0); - _nodeMap[std::pair<int,int>(0,0)]=_root; listPrev.push_back(_root); } else @@ -699,7 +684,6 @@ void Tree::createMixedTree(int maxDepth, int numChildMax, int leafPerNode, int n { np->childs[j] = new TreeNode(depth,inode,np,j); } - _nodeMap[std::pair<int,int>(depth,inode)]=np->childs[j]; listCur.push_back(np->childs[j]); inode++; } @@ -729,7 +713,6 @@ void Tree::createPerfectTree(int maxDepth, int numChildRegular, int numChildLeaf if (depth == 0) { _root = new TreeNode(0,0,NULL,0); - _nodeMap[std::pair<int,int>(0,0)]=_root; listPrev.push_back(_root); } else @@ -756,7 +739,6 @@ void Tree::createPerfectTree(int maxDepth, int numChildRegular, int numChildLeaf { np->childs[j] = new TreeNode(depth,inode,np,j); } - _nodeMap[std::pair<int,int>(depth,inode)]=np->childs[j]; listCur.push_back(np->childs[j]); inode++; } @@ -787,7 +769,6 @@ void Tree::createPerfectTree(int maxDepth, int numChilds, int numDirVars,const s if (depth == 0) { _root = new TreeNode(0,0,NULL,0); - _nodeMap[std::pair<int,int>(0,0)]=_root; listPrev.push_back(_root); } else @@ -809,7 +790,6 @@ void Tree::createPerfectTree(int maxDepth, int numChilds, int numDirVars,const s { np->childs[j] = new TreeNode(depth,inode,np,j); } - _nodeMap[std::pair<int,int>(depth,inode)]=np->childs[j]; listCur.push_back(np->childs[j]); inode++; } @@ -865,7 +845,7 @@ void Tree::loadDataFromFile(std::string filename) { // clear tree clear(); - _nodeMap.clear(); + std::map<std::pair<int,int>, TreeNode* > nodeMap; printf("loading data from file %s\n",filename.c_str()); // FILE* fp = fopen(filename.c_str(),"r"); @@ -889,20 +869,20 @@ void Tree::loadDataFromFile(std::string filename) std::vector<int> allChilds(numChild); for (int i=0; i< numChild; i++) { - fscanf(fp, "%d",&allChilds[i]); + int fl=fscanf(fp, "%d",&allChilds[i]); printf("child %d location %d\n",i,allChilds[i]); } double weight; - fscanf(fp, "%lf",&weight); + int fl=fscanf(fp, "%lf",&weight); printf("weight = %f\n",weight); std::vector<double> direction(numDir); for (int i=0; i< numDir; i++) { - fscanf(fp, "%lf",&direction[i]); + fl=fscanf(fp, "%lf",&direction[i]); printf("direction %d value %f\n",i,direction[i]); } char what[256]; - fscanf(fp, "%s", what); + fl=fscanf(fp, "%s", what); printf("aftype = %s\n",what); // create node @@ -926,7 +906,7 @@ void Tree::loadDataFromFile(std::string filename) node->childs.clear(); } _root = node; - _nodeMap[std::pair<int,int>(depth,location)]=node; + nodeMap[std::pair<int,int>(depth,location)]=node; } else { @@ -937,7 +917,7 @@ void Tree::loadDataFromFile(std::string filename) node->depth = depth; node->location = location; node->childOrder = childOrder; - node->parent = _nodeMap[std::pair<int,int>(depth-1,parentIdLoc)]; + node->parent = nodeMap[std::pair<int,int>(depth-1,parentIdLoc)]; node->parent->childs[childOrder] = node; node->matIndex = matIndex; if (numChild > 0) @@ -950,7 +930,7 @@ void Tree::loadDataFromFile(std::string filename) } printf("node %d %d:\n",depth,location); node->parent->printData(); - _nodeMap[std::pair<int,int>(depth,location)]=node; + nodeMap[std::pair<int,int>(depth,location)]=node; } } fclose(fp); @@ -1027,7 +1007,7 @@ bool Tree::removeZeroLeaves(double tol, bool iteration) otherChild->parent = parent->parent; otherChild->childOrder = parent->childOrder; grandParent->childs[parent->childOrder] = otherChild; - Msg::Info("otherChild weight = %e parent weight = %e",otherChild->weight,parent->weight); + Msg::Info("child weight = %e otherChild weight = %e parent weight = %e",leaf->af->getVal(leaf->weight),otherChild->af->getVal(otherChild->weight),parent->weight); } else if (parent->childs.size() > 2) { @@ -1263,6 +1243,15 @@ bool Tree::treeMerging() } } }; + + Msg::Info("\nSTART merging "); + parent->printData(); + printf("merge loc (number of phases = %ld): %d (mat %d) ",parent->childs.size(),sameMat[0].second, (sameMat[0].first)->matIndex); + for (std::set<int>::iterator itpos = position.begin(); itpos != position.end(); itpos++) + { + printf(" %d (mat %d)",*itpos,parent->childs[*itpos]->matIndex); + } + printf("\n----------------------\n"); parent->childs = newChilds; for (int j=0; j <parent->childs.size(); j++) { @@ -1270,6 +1259,7 @@ bool Tree::treeMerging() } parent->direction = newNormals; parent->printData(); + Msg::Info("DONE merging\n"); // break; } @@ -1281,6 +1271,110 @@ bool Tree::treeMerging() return false; }; +bool Tree::removeZeroBranches(double tol) +{ + int depth = 0; + bool ok = false; + while (true) + { + std::map<int,std::vector<TreeNode*> > allNodesByDepth; + getAllNodesByDepth(allNodesByDepth); + std::map<int,std::vector<TreeNode*> >::iterator itF= allNodesByDepth.find(depth); + if ( itF!= allNodesByDepth.end()) + { + std::vector<TreeNode*>& depthNodes = itF->second; + for (int i=0; i<depthNodes.size(); i++) + { + TreeNode* node = depthNodes[i]; + double Wt = getNonNegativeWeight(node); + // + bool removed = false; + std::vector<TreeNode*> newChilds; + std::vector<double> newDirection; + std::vector<int> removeLoc; + // + int numberNormal = node->childs.size()-1; + int totalNumberDirVars = node->direction.size(); + int numberVarsPerNormal = totalNumberDirVars/numberNormal; + + for (int j=0; j< node->childs.size(); j++) + { + double Wchild = getNonNegativeWeight(node->childs[j]); + if (Wchild < tol* Wt) + { + removed = true; + removeLoc.push_back(j); + } + else + { + newChilds.push_back(node->childs[j]); + if (newChilds.size() > 1) + { + if (numberVarsPerNormal == 1) + { + newDirection.push_back(node->direction[j-1]); + } + else if (numberVarsPerNormal == 2) + { + newDirection.push_back(node->direction[(j-1)*numberVarsPerNormal]); + newDirection.push_back(node->direction[(j-1)*numberVarsPerNormal+1]); + } + else + { + Msg::Error("numberVarsPerNormal = %d is not correct !!!",numberVarsPerNormal); + Msg::Exit(0); + } + } + } + } + if (removed) + { + Msg::Info("\nSTART removing "); + node->printData(); + ok = true; + printf("pos (number of phases= %ld): ",node->childs.size()); + for (int j=0; j< removeLoc.size(); j++) + { + printf("%d (f=%e) ",removeLoc[j],getNonNegativeWeight(node->childs[removeLoc[j]])/Wt); + } + printf("\n"); + if (newChilds.size() >1) + { + node->childs = newChilds; + for (int j=0; j <node->childs.size(); j++) + { + node->childs[j]->childOrder = j; + } + node->direction = newDirection; + Msg::Info("------------------"); + node->printData(); + } + else + { + TreeNode* grandParent = node->parent; + TreeNode* child = newChilds[0]; + // + child->depth = node->depth; + child->location = node->location; + child->parent = node->parent; + child->childOrder = node->childOrder; + grandParent->childs[node->childOrder] = child; + Msg::Info("------------------"); + Msg::Info("node is replaced by child"); + } + Msg::Info("DONE removing\n"); + } + } + } + else + { + break; + } + depth++; + }; + return ok; +}; + bool Tree::removeLeavesWithZeroContribution(double tol) { int iterRemove = 0; @@ -1294,7 +1388,7 @@ bool Tree::removeLeavesWithZeroContribution(double tol) Msg::Error("removeLeavesWithZeroContribution: maximal number of iterations (%d) reaches", maxIterRemove); Msg::Exit(0); } - bool okLeaf = removeZeroLeaves(tol); + bool okLeaf = removeZeroBranches(tol); bool okMerge = treeMerging(); // check if (okLeaf || okLeaf) @@ -1396,7 +1490,7 @@ void Tree::printTree() const nodeContainer allLeaves; getAllLeaves(allLeaves); - printf("all leaves = %d\n",allLeaves.size()); + printf("all leaves = %ld\n",allLeaves.size()); for (int i=0; i< allLeaves.size(); i++) { allLeaves[i]->printData(); @@ -1412,7 +1506,7 @@ void Tree::printTreeInteraction(const std::string fname, bool colorMat, int dir) listPrev.push_back(_root); nodeContainer allLeaves; getAssociatedLeavesForNode(_root,allLeaves); - printf("*|%d|*\n",allLeaves.size()); + printf("*|%ld|*\n",allLeaves.size()); while (true) { listCur.clear(); @@ -1429,22 +1523,22 @@ void Tree::printTreeInteraction(const std::string fname, bool colorMat, int dir) { if (j==0) { - printf("%d(%d)",allLeaves.size(),nj->matIndex); + printf("%ld(%d)",allLeaves.size(),nj->matIndex); } else { - printf("-%d(%d)",allLeaves.size(),nj->matIndex); + printf("-%ld(%d)",allLeaves.size(),nj->matIndex); } } else { if (j==0) { - printf("%d",allLeaves.size()); + printf("%ld",allLeaves.size()); } else { - printf("-%d",allLeaves.size()); + printf("-%ld",allLeaves.size()); } listCur.push_back(nj); } @@ -1596,11 +1690,11 @@ void Tree::getRefToAllLeaves(std::vector<TreeNode*>& allLeaves) }; }; -void Tree::getAllNodesByDepth(std::map<int,nodeContainer >& allNodes) const +void Tree::getAllNodesByDepth(std::map<int,std::vector<TreeNode*> >& allNodes) { allNodes.clear(); if (_root == NULL) return; - nodeContainer listPrev, listCur; + std::vector<TreeNode*> listPrev, listCur; listPrev.push_back(_root); // allNodes[_root->depth].push_back(_root); @@ -1609,7 +1703,7 @@ void Tree::getAllNodesByDepth(std::map<int,nodeContainer >& allNodes) const listCur.clear(); for (int i =0; i< listPrev.size(); i++) { - const TreeNode* np = listPrev[i]; + TreeNode* np = listPrev[i]; for (int j=0; j< np->childs.size(); j++) { listCur.push_back(np->childs[j]); @@ -1624,24 +1718,42 @@ void Tree::getAllNodesByDepth(std::map<int,nodeContainer >& allNodes) const }; }; -const TreeNode* Tree::getNode(int depth, int pos) const +void Tree::getAllNodesByDepth(std::map<int,nodeContainer >& allNodes) const { - std::map<std::pair<int,int>, TreeNode* >::const_iterator it = _nodeMap.find(std::pair<int,int>(depth,pos)); - if (it == _nodeMap.end()) + allNodes.clear(); + if (_root == NULL) return; + nodeContainer listPrev, listCur; + listPrev.push_back(_root); + // + allNodes[_root->depth].push_back(_root); + while (true) { - printf("node in depth %d pos %d is not found \n",depth,pos); - return NULL; - } - return it->second; + listCur.clear(); + for (int i =0; i< listPrev.size(); i++) + { + const TreeNode* np = listPrev[i]; + for (int j=0; j< np->childs.size(); j++) + { + listCur.push_back(np->childs[j]); + allNodes[np->childs[j]->depth].push_back(np->childs[j]); + } + } + listPrev = listCur; + if (listCur.size() == 0) + { + break; + } + }; }; + void Tree::printBackTrackingPath(const TreeNode* node) const { nodeContainer allNodes; getBackTrackingPath(node,allNodes); printf("-----------------------------\n"); node->printData(); - printf("backtracking laves to node: %d \n",allNodes.size()); + printf("backtracking laves to node: %ld \n",allNodes.size()); for (int i=0; i< allNodes.size(); i++) { allNodes[i]->printData(); @@ -1656,7 +1768,7 @@ void Tree::printAssociatedLeavesForNode(const TreeNode* node) const printf("-----------------------------\n"); node->printData(); - printf("associated laves to node: %d \n",allLeaves.size()); + printf("associated laves to node: %ld \n",allLeaves.size()); for (int i=0; i< allLeaves.size(); i++) { allLeaves[i]->printData(); @@ -1676,6 +1788,54 @@ void Tree::getBackTrackingPath(const TreeNode* node, nodeContainer& allNodes) co }; }; +void Tree::getAssociatedLeavesForNode(TreeNode* node, std::vector<TreeNode*>& allLeaves) +{ + allLeaves.clear(); + std::vector<TreeNode*> listPrev, listCur; + listPrev.push_back(node); + while (true) + { + listCur.clear(); + for (int i =0; i< listPrev.size(); i++) + { + TreeNode* np = listPrev[i]; + if (np->childs.size() == 0) + { + allLeaves.push_back(np); + } + for (int j=0; j< np->childs.size(); j++) + { + listCur.push_back(np->childs[j]); + } + } + listPrev = listCur; + if (listCur.size() == 0) + { + break; + } + }; +}; + +double Tree::getNonNegativeWeight(const TreeNode* node) const +{ + if (node->childs.size() == 0) + { + return node->af->getVal(node->weight); + } + else + { + nodeContainer allLeaves; + getAssociatedLeavesForNode(node,allLeaves); + double W = 0; + for (int i=0; i< allLeaves.size(); i++) + { + const TreeNode* leaf = allLeaves[i]; + W += leaf->af->getVal(leaf->weight); + }; + return W; + } +} + void Tree::getAssociatedLeavesForNode(const TreeNode* node, nodeContainer& allLeaves) const { allLeaves.clear(); @@ -1763,7 +1923,7 @@ void Tree::printAllRegularNodes() const nodeContainer allNodes; getAllRegularNodes(allNodes); printf("-----------------------------\n"); - printf("all regular nodes: %d \n",allNodes.size()); + printf("all regular nodes: %ld \n",allNodes.size()); for (int i=0; i< allNodes.size(); i++) { allNodes[i]->printData(); @@ -1926,7 +2086,7 @@ void Tree::initialize(bool rand) // std::vector<TreeNode*> allNodes; getRefToAllNodes(allNodes); - printf("num of nodes = %d\n",allNodes.size()); + printf("num of nodes = %ld\n",allNodes.size()); for (int i=0; i< allNodes.size(); i++) { TreeNode& node = *(allNodes[i]); diff --git a/NonLinearSolver/modelReduction/Tree.h b/NonLinearSolver/modelReduction/Tree.h index e5d7abd6622e4ddd268951038980af00e6afd7ac..b2a94811333a323b9c9cd4bef6896bf3af0ee6f3 100644 --- a/NonLinearSolver/modelReduction/Tree.h +++ b/NonLinearSolver/modelReduction/Tree.h @@ -61,7 +61,6 @@ class Tree void clear(); void initialize(bool rand=true); const TreeNode* getRootNode() const {return _root;}; - const TreeNode* getNode(int depth, int pos) const; void assignMaterialLaws(int numPhases); @@ -78,6 +77,7 @@ class Tree void printNodeFraction() const; bool removeLeavesWithZeroContribution(double tol=1e-6); + bool removeZeroBranches(double tol); bool removeZeroLeaves(double tol, bool iteration=true); bool treeMerging(); @@ -86,9 +86,11 @@ class Tree void getAllLeaves(nodeContainer& allLeaves) const; void getAllNodes(nodeContainer& allNodes) const; void getAllNodesByDepth(std::map<int,nodeContainer>& allNodes) const; + void getAllNodesByDepth(std::map<int,std::vector<TreeNode*> >& allNodes); void getAssociatedLeavesForNode(const TreeNode* node, nodeContainer& allLeaves) const; + void getAssociatedLeavesForNode(TreeNode* node, std::vector<TreeNode*>& allLeaves); void getBackTrackingPath(const TreeNode* node, nodeContainer& path) const; - + double getNonNegativeWeight(const TreeNode* node) const; void getRefToAllLeaves(std::vector<TreeNode*>& allLeaves); void getRefToAllNodes(std::vector<TreeNode*>& allNodes); #endif //SWIG @@ -96,7 +98,6 @@ class Tree protected: #ifndef SWIG TreeNode* _root; - std::map<std::pair<int,int>, TreeNode* > _nodeMap; #endif //SWIG }; #endif //_TREE_H_