1 /*---------------------------------------------------------------------------*\
3 \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
5 \\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
7 -------------------------------------------------------------------------------
9 This file is part of OpenFOAM.
11 OpenFOAM is free software; you can redistribute it and/or modify it
12 under the terms of the GNU General Public License as published by the
13 Free Software Foundation; either version 2 of the License, or (at your
14 option) any later version.
16 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
21 You should have received a copy of the GNU General Public License
22 along with OpenFOAM; if not, write to the Free Software Foundation,
23 Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
25 \*---------------------------------------------------------------------------*/
27 #include "slidingInterface.H"
28 #include "polyTopoChanger.H"
30 #include "primitiveMesh.H"
31 #include "polyTopoChange.H"
32 #include "addToRunTimeSelectionTable.H"
33 #include "triPointRef.H"
36 // Index of debug signs:
37 // p - adjusting a projection point
38 // * - adjusting edge intersection
40 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
44 defineTypeNameAndDebug(slidingInterface, 0);
45 addToRunTimeSelectionTable
55 const char* Foam::NamedEnum<Foam::slidingInterface::typeOfMatch, 2>::names[] =
62 const Foam::NamedEnum<Foam::slidingInterface::typeOfMatch, 2>
63 Foam::slidingInterface::typeOfMatchNames_;
66 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
68 void Foam::slidingInterface::checkDefinition()
70 const polyMesh& mesh = topoChanger().mesh();
74 !masterFaceZoneID_.active()
75 || !slaveFaceZoneID_.active()
76 || !cutPointZoneID_.active()
77 || !cutFaceZoneID_.active()
78 || !masterPatchID_.active()
79 || !slavePatchID_.active()
84 "void slidingInterface::checkDefinition()"
85 ) << "Not all zones and patches needed in the definition "
86 << "have been found. Please check your mesh definition."
90 // Check the sizes and set up state
93 mesh.faceZones()[masterFaceZoneID_.index()].size() == 0
94 || mesh.faceZones()[slaveFaceZoneID_.index()].size() == 0
97 FatalErrorIn("void slidingInterface::checkDefinition()")
98 << "Master or slave face zone contain no faces. "
99 << "Please check your mesh definition."
100 << abort(FatalError);
105 Pout<< "Sliding interface object " << name() << " :" << nl
106 << " master face zone: " << masterFaceZoneID_.index() << nl
107 << " slave face zone: " << slaveFaceZoneID_.index() << endl;
112 void Foam::slidingInterface::clearOut() const
114 clearPointProjection();
115 clearAttachedAddressing();
120 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
123 // Construct from components
124 Foam::slidingInterface::slidingInterface
128 const polyTopoChanger& mme,
129 const word& masterFaceZoneName,
130 const word& slaveFaceZoneName,
131 const word& cutPointZoneName,
132 const word& cutFaceZoneName,
133 const word& masterPatchName,
134 const word& slavePatchName,
135 const typeOfMatch tom,
136 const bool coupleDecouple,
137 const intersection::algorithm algo
140 polyMeshModifier(name, index, mme, true),
144 mme.mesh().faceZones()
149 mme.mesh().faceZones()
154 mme.mesh().pointZones()
159 mme.mesh().faceZones()
164 mme.mesh().boundaryMesh()
169 mme.mesh().boundaryMesh()
172 coupleDecouple_(coupleDecouple),
174 projectionAlgo_(algo),
176 cutFaceMasterPtr_(NULL),
177 cutFaceSlavePtr_(NULL),
178 masterFaceCellsPtr_(NULL),
179 slaveFaceCellsPtr_(NULL),
180 masterStickOutFacesPtr_(NULL),
181 slaveStickOutFacesPtr_(NULL),
182 retiredPointMapPtr_(NULL),
183 cutPointEdgePairMapPtr_(NULL),
184 slavePointPointHitsPtr_(NULL),
185 slavePointEdgeHitsPtr_(NULL),
186 slavePointFaceHitsPtr_(NULL),
187 masterPointEdgeHitsPtr_(NULL),
188 projectedSlavePointsPtr_(NULL)
196 "Foam::slidingInterface::slidingInterface\n"
198 " const word& name,\n"
199 " const label index,\n"
200 " const polyTopoChanger& mme,\n"
201 " const word& masterFaceZoneName,\n"
202 " const word& slaveFaceZoneName,\n"
203 " const word& cutFaceZoneName,\n"
204 " const word& cutPointZoneName,\n"
205 " const word& masterPatchName,\n"
206 " const word& slavePatchName,\n"
207 " const typeOfMatch tom,\n"
208 " const bool coupleDecouple\n"
210 ) << "Creation of a sliding interface from components "
211 << "in attached state not supported."
212 << abort(FatalError);
216 calcAttachedAddressing();
221 // Construct from components
222 Foam::slidingInterface::slidingInterface
225 const dictionary& dict,
227 const polyTopoChanger& mme
230 polyMeshModifier(name, index, mme, Switch(dict.lookup("active"))),
233 dict.lookup("masterFaceZoneName"),
234 mme.mesh().faceZones()
238 dict.lookup("slaveFaceZoneName"),
239 mme.mesh().faceZones()
243 dict.lookup("cutPointZoneName"),
244 mme.mesh().pointZones()
248 dict.lookup("cutFaceZoneName"),
249 mme.mesh().faceZones()
253 dict.lookup("masterPatchName"),
254 mme.mesh().boundaryMesh()
258 dict.lookup("slavePatchName"),
259 mme.mesh().boundaryMesh()
261 matchType_(typeOfMatchNames_.read((dict.lookup("typeOfMatch")))),
262 coupleDecouple_(dict.lookup("coupleDecouple")),
263 attached_(dict.lookup("attached")),
266 intersection::algorithmNames_.read(dict.lookup("projection"))
269 cutFaceMasterPtr_(NULL),
270 cutFaceSlavePtr_(NULL),
271 masterFaceCellsPtr_(NULL),
272 slaveFaceCellsPtr_(NULL),
273 masterStickOutFacesPtr_(NULL),
274 slaveStickOutFacesPtr_(NULL),
275 retiredPointMapPtr_(NULL),
276 cutPointEdgePairMapPtr_(NULL),
277 slavePointPointHitsPtr_(NULL),
278 slavePointEdgeHitsPtr_(NULL),
279 slavePointFaceHitsPtr_(NULL),
280 masterPointEdgeHitsPtr_(NULL),
281 projectedSlavePointsPtr_(NULL)
285 // If the interface is attached, the master and slave face zone addressing
286 // needs to be read in; otherwise it will be created
291 Pout<< "slidingInterface::slidingInterface(...) "
292 << " for object " << name << " : "
293 << "Interface attached. Reading master and slave face zones "
294 << "and retired point lookup." << endl;
297 // The face zone addressing is written out in the definition dictionary
298 masterFaceCellsPtr_ = new labelList(dict.lookup("masterFaceCells"));
299 slaveFaceCellsPtr_ = new labelList(dict.lookup("slaveFaceCells"));
301 masterStickOutFacesPtr_ =
302 new labelList(dict.lookup("masterStickOutFaces"));
303 slaveStickOutFacesPtr_ =
304 new labelList(dict.lookup("slaveStickOutFaces"));
306 retiredPointMapPtr_ = new Map<label>(dict.lookup("retiredPointMap"));
307 cutPointEdgePairMapPtr_ =
308 new Map<Pair<edge> >(dict.lookup("cutPointEdgePairMap"));
312 calcAttachedAddressing();
317 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
319 Foam::slidingInterface::~slidingInterface()
325 void Foam::slidingInterface::clearAddressing() const
327 deleteDemandDrivenData(cutFaceMasterPtr_);
328 deleteDemandDrivenData(cutFaceSlavePtr_);
332 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
334 const Foam::faceZoneID& Foam::slidingInterface::masterFaceZoneID() const
336 return masterFaceZoneID_;
340 const Foam::faceZoneID& Foam::slidingInterface::slaveFaceZoneID() const
342 return slaveFaceZoneID_;
346 bool Foam::slidingInterface::changeTopology() const
350 // Always changes. If not attached, project points
353 Pout<< "bool slidingInterface::changeTopology() const "
354 << "for object " << name() << " : "
355 << "Couple-decouple mode." << endl;
372 && !topoChanger().mesh().changing()
375 // If the mesh is not moving or morphing and the interface is
376 // already attached, the topology will not change
381 // Check if the motion changes point projection
382 return projectPoints();
387 void Foam::slidingInterface::setRefinement(polyTopoChange& ref) const
393 // Attached, coupling
394 decoupleInterface(ref);
398 // Detached, coupling
399 coupleInterface(ref);
409 // Clear old coupling data
413 coupleInterface(ref);
420 void Foam::slidingInterface::modifyMotionPoints(pointField& motionPoints) const
424 Pout<< "void slidingInterface::modifyMotionPoints("
425 << "pointField& motionPoints) const for object " << name() << " : "
426 << "Adjusting motion points." << endl;
429 const polyMesh& mesh = topoChanger().mesh();
431 // Get point from the cut zone
432 const labelList& cutPoints = mesh.pointZones()[cutPointZoneID_.index()];
434 if (cutPoints.size() > 0 && !projectedSlavePointsPtr_)
440 const pointField& projectedSlavePoints = *projectedSlavePointsPtr_;
442 const Map<label>& rpm = retiredPointMap();
444 const Map<Pair<edge> >& cpepm = cutPointEdgePairMap();
446 const Map<label>& slaveZonePointMap =
447 mesh.faceZones()[slaveFaceZoneID_.index()]().meshPointMap();
449 const primitiveFacePatch& masterPatch =
450 mesh.faceZones()[masterFaceZoneID_.index()]();
451 const edgeList& masterEdges = masterPatch.edges();
452 const pointField& masterLocalPoints = masterPatch.localPoints();
454 const primitiveFacePatch& slavePatch =
455 mesh.faceZones()[slaveFaceZoneID_.index()]();
456 const edgeList& slaveEdges = slavePatch.edges();
457 const pointField& slaveLocalPoints = slavePatch.localPoints();
458 const vectorField& slavePointNormals = slavePatch.pointNormals();
460 forAll (cutPoints, pointI)
462 // Try to find the cut point in retired points
463 Map<label>::const_iterator rpmIter = rpm.find(cutPoints[pointI]);
465 if (rpmIter != rpm.end())
472 // Cut point is a retired point
473 motionPoints[cutPoints[pointI]] =
474 projectedSlavePoints[slaveZonePointMap.find(rpmIter())()];
478 // A cut point is not a projected slave point. Therefore, it
479 // must be an edge-to-edge intersection.
481 Map<Pair<edge> >::const_iterator cpepmIter =
482 cpepm.find(cutPoints[pointI]);
484 if (cpepmIter != cpepm.end())
486 // Pout << "Need to re-create hit for point " << cutPoints[pointI] << " lookup: " << cpepmIter() << endl;
489 // The edge cutting code is repeated in
490 // slidingInterface::coupleInterface. This is done for
491 // efficiency reasons and avoids multiple creation of
492 // cutting planes. Please update both simultaneously.
494 const edge& globalMasterEdge = cpepmIter().first();
496 const label curMasterEdgeIndex =
497 masterPatch.whichEdge
501 masterPatch.whichPoint
503 globalMasterEdge.start()
505 masterPatch.whichPoint
507 globalMasterEdge.end()
512 const edge& cme = masterEdges[curMasterEdgeIndex];
513 // Pout << "curMasterEdgeIndex: " << curMasterEdgeIndex << " cme: " << cme << endl;
514 const edge& globalSlaveEdge = cpepmIter().second();
516 const label curSlaveEdgeIndex =
521 slavePatch.whichPoint
523 globalSlaveEdge.start()
525 slavePatch.whichPoint
527 globalSlaveEdge.end()
532 const edge& curSlaveEdge = slaveEdges[curSlaveEdgeIndex];
533 // Pout << "curSlaveEdgeIndex: " << curSlaveEdgeIndex << " curSlaveEdge: " << curSlaveEdge << endl;
534 const point& a = projectedSlavePoints[curSlaveEdge.start()];
535 const point& b = projectedSlavePoints[curSlaveEdge.end()];
540 slaveLocalPoints[curSlaveEdge.start()]
541 + slavePointNormals[curSlaveEdge.start()]
542 + slaveLocalPoints[curSlaveEdge.end()]
543 + slavePointNormals[curSlaveEdge.end()]
547 plane cutPlane(a, b, c);
549 linePointRef curSlaveLine =
550 curSlaveEdge.line(slaveLocalPoints);
551 const scalar curSlaveLineMag = curSlaveLine.mag();
554 cutPlane.lineIntersect
556 cme.line(masterLocalPoints)
561 cutOnMaster > edgeEndCutoffTol_
562 && cutOnMaster < 1.0 - edgeEndCutoffTol_
565 // Master is cut, check the slave
566 point masterCutPoint =
567 masterLocalPoints[cme.start()]
568 + cutOnMaster*cme.vec(masterLocalPoints);
571 curSlaveLine.nearestDist(masterCutPoint);
575 // Strict checking of slave cut to avoid capturing
581 - curSlaveLine.start()
582 ) & curSlaveLine.vec()
583 )/sqr(curSlaveLineMag);
585 // Calculate merge tolerance from the
586 // target edge length
588 edgeCoPlanarTol_*mag(b - a);
592 cutOnSlave > edgeEndCutoffTol_
593 && cutOnSlave < 1.0 - edgeEndCutoffTol_
594 && slaveCut.distance() < mergeTol
597 // Cut both master and slave.
598 motionPoints[cutPoints[pointI]] =
604 Pout<< "Missed slave edge!!! This is an error. "
606 << cme.line(masterLocalPoints)
607 << " slave edge: " << curSlaveLine
608 << " point: " << masterCutPoint
613 - curSlaveLine.start()
614 ) & curSlaveLine.vec()
615 )/sqr(curSlaveLineMag)
621 Pout<< "Missed master edge!!! This is an error"
629 "void slidingInterface::modifyMotionPoints"
630 "(pointField&) const"
631 ) << "Cut point " << cutPoints[pointI]
632 << " not recognised as either the projected "
633 << "or as intersection point. Error in point "
634 << "projection or data mapping"
635 << abort(FatalError);
647 void Foam::slidingInterface::updateMesh(const mapPolyMesh& m)
651 Pout<< "void slidingInterface::updateMesh(const mapPolyMesh& m)"
652 << " const for object " << name() << " : "
653 << "Updating topology." << endl;
656 // Mesh has changed topologically. Update local topological data
657 const polyMesh& mesh = topoChanger().mesh();
659 masterFaceZoneID_.update(mesh.faceZones());
660 slaveFaceZoneID_.update(mesh.faceZones());
661 cutPointZoneID_.update(mesh.pointZones());
662 cutFaceZoneID_.update(mesh.faceZones());
664 masterPatchID_.update(mesh.boundaryMesh());
665 slavePatchID_.update(mesh.boundaryMesh());
668 //Pout<< "**Mj: disabled updating." << endl;
672 // calcAttachedAddressing();
676 // renumberAttachedAddressing(m);
681 const Foam::pointField& Foam::slidingInterface::pointProjection() const
683 if (!projectedSlavePointsPtr_)
688 return *projectedSlavePointsPtr_;
692 void Foam::slidingInterface::write(Ostream& os) const
694 os << nl << type() << nl
696 << masterFaceZoneID_.name() << nl
697 << slaveFaceZoneID_.name() << nl
698 << cutPointZoneID_.name() << nl
699 << cutFaceZoneID_.name() << nl
700 << masterPatchID_.name() << nl
701 << slavePatchID_.name() << nl
702 << typeOfMatchNames_[matchType_] << nl
703 << coupleDecouple_ << nl
704 << attached_ << endl;
708 void Foam::slidingInterface::writeDict(Ostream& os) const
710 os << nl << name() << nl << token::BEGIN_BLOCK << nl
711 << " type " << type() << token::END_STATEMENT << nl
712 << " masterFaceZoneName " << masterFaceZoneID_.name()
713 << token::END_STATEMENT << nl
714 << " slaveFaceZoneName " << slaveFaceZoneID_.name()
715 << token::END_STATEMENT << nl
716 << " cutPointZoneName " << cutPointZoneID_.name()
717 << token::END_STATEMENT << nl
718 << " cutFaceZoneName " << cutFaceZoneID_.name()
719 << token::END_STATEMENT << nl
720 << " masterPatchName " << masterPatchID_.name()
721 << token::END_STATEMENT << nl
722 << " slavePatchName " << slavePatchID_.name()
723 << token::END_STATEMENT << nl
724 << " typeOfMatch " << typeOfMatchNames_[matchType_]
725 << token::END_STATEMENT << nl
726 << " coupleDecouple " << coupleDecouple_
727 << token::END_STATEMENT << nl
728 << " projection " << intersection::algorithmNames_[projectionAlgo_]
729 << token::END_STATEMENT << nl
730 << " attached " << attached_
731 << token::END_STATEMENT << nl
732 << " active " << active()
733 << token::END_STATEMENT << nl;
737 masterFaceCellsPtr_->writeEntry("masterFaceCells", os);
738 slaveFaceCellsPtr_->writeEntry("slaveFaceCells", os);
739 masterStickOutFacesPtr_->writeEntry("masterStickOutFaces", os);
740 slaveStickOutFacesPtr_->writeEntry("slaveStickOutFaces", os);
742 os << " retiredPointMap " << retiredPointMap()
743 << token::END_STATEMENT << nl
744 << " cutPointEdgePairMap " << cutPointEdgePairMap()
745 << token::END_STATEMENT << nl;
748 os << token::END_BLOCK << endl;
752 // ************************************************************************* //