Merge pull request #380 from ONSAS/jorge

vtk quad hexa
ONSAS · Jan 20, 2022 · 59da993 · 59da993
2 parents c32f32f + 72f7d57
commit 59da993
Show file tree

Hide file tree

Showing 9 changed files with 114 additions and 143 deletions.
diff --git a/README.md b/README.md
@@ -15,7 +15,6 @@
 1. [Contact](#contact)
 
 ## About ONSAS <a name="aboutonsas"></a>
-------
 
 ### What is ONSAS?
 
@@ -57,7 +56,6 @@ The user should follow these steps to install and run onsas:
 1. Open GNU-Octave, move to the _examples_ directory and run one of the examples.
 
 ## Contributions <a name="contributions"></a>
-------
 
 ### Authors of code
 The following authors collaborated in various tasks including: design, development and testing of the code, with specific contributions given by the commits history: [**Jorge M. Pérez Zerpa**](https://www.fing.edu.uy/~jorgepz) <sup>1</sup>, **J. Bruno Bazzano**<sup>1,2</sup>, [**Joaquín Viera**](https://www.researchgate.net/profile/Joaquin_Viera_Sosa) <sup>1</sup>, **Mauricio Vanzulli** <sup>3</sup> and [**Marcelo Forets**](https://scholar.google.fr/citations?user=XSJzDEsAAAAJ&hl=en)<sup>4</sup>.
@@ -78,6 +76,5 @@ Affiliations:
 The complete list of contributions and acknowledgments is available at the documentation site.
 
 ## Contact <a name="contact"></a>
-------
 
-If you have any question you can send an e-mail to _jorgepz[AT]fing.edu.uy_ or [![Join the chat at https://gitter.im/onsas_/community](https://badges.gitter.im/onsas_/community.svg)](https://gitter.im/onsas_/community?utm_source=badge&utm_medium=badge&utm_campaign=pr-badge&utm_content=badge) .
+If you have any question you can send an e-mail to _jorgepz [AT] fing.edu.uy_ or [![Join the chat at https://gitter.im/onsas_/community](https://badges.gitter.im/onsas_/community.svg)](https://gitter.im/onsas_/community?utm_source=badge&utm_medium=badge&utm_campaign=pr-badge&utm_content=badge) .
diff --git a/about_docs.md → docs/about_docs.md b/about_docs.md → docs/about_docs.md
diff --git a/examples/beamTrussPendulum/onsasExample_beamTrussPendulum.m b/examples/beamTrussPendulum/onsasExample_beamTrussPendulum.m
@@ -3,12 +3,17 @@
 %mdProblem name:
 otherParams.problemName = 'beamTrussPendulum' ;
 addpath( genpath( [ pwd '/../../src'] ) );
-%mdTuss element geometrical properties :
-Et = 10e11 ; nut = 0.3 ;  rhot = 65.6965 ;
-At = .1 ; dt = sqrt(4*At/pi);   lt = 3.0443 ;  
+%mdTuss element material and geometrical scalar properties
+Et = 200e9 ; nut = 0.3 ;  rhot = 8e3 ;
+dt = .02 ; At = pi*dt^2/4 ; lt = 1 ;
+%
 %mdFrame element geometrical properties :
-Ef = Et/300000*7; nuf = 0.3;  rhof = rhot;
-df = dt; Ab = pi*df^2/4;  lf = lt; If = pi*df^4/64 ;  
+Ef = Et/20 ; nuf = 0.3;  rhof = rhot;
+df = dt*2; Ab = pi*df^2/4;  lf = lt; If = pi*df^4/64 ;
+%
+%mdScalar parameters of the mesh:
+numElemF = 10; numElemT = 1;
+%
 %md### MEBI parameters
 %md
 %md### materials
@@ -28,7 +33,7 @@
 %md
 %mdTruss:
 elements(2).elemType = 'truss'          ;
-elements(2).elemTypeGeometry = [ 3 dt ] ;
+elements(2).elemTypeGeometry = [ 2 dt dt ] ;
 elements(2).elemTypeParams = 1          ;
 %md
 %mdFrame:
@@ -51,42 +56,36 @@
 initialConds                = struct() ;
 %md### mesh parameters
 %md
-%mdScalar parameters of the mesh:
-numElemF = 10; numElemT = 1;
-%md
 %mdnodesCoords:
-mesh.nodesCoords = [  ( 0:(numElemF)  )' * lf/numElemF    zeros(numElemF + 1,1)   zeros(numElemF + 1,1);	 
-		                  ( lf + (1:numElemT)' * lt/numElemT )  zeros(numElemT,1)       zeros(numElemT,1)    ];
+mesh.nodesCoords = [  ( 0:(lf/numElemF):lf)'        zeros(numElemF + 1, 1)  zeros(numElemF + 1,1)  ; ...
+		                  ( lf+(1:(lt/numElemT):lt)' )  zeros(numElemT    , 1)  zeros(numElemT,1)      ] ;
 %md
 %mdConec matrix:
 %mdnodes conectivity:
 mesh.conecCell = { } ;
 mesh.conecCell{ 1, 1 } = [ 0 1 1 0  1                   ] ;
-mesh.conecCell{ 2, 1 } = [ 0 1 2 0  numElemF + 1        ] ;
-mesh.conecCell{ 3, 1 } = [ 0 1 2 0  numElemF + numElemT ] ;
+mesh.conecCell{ 2, 1 } = [ 0 1 2 0  numElemF+1+numElemT ] ;
+%
 %mdelements conectivity:
 auxConecElem  = [ %MEBI frame elements
-                  [ (ones(numElemF,1)*2 )    (ones(numElemF,1)*3)      (zeros(numElemF,1))    (zeros(numElemF,1)) ... 
+                  [ (ones(numElemF,1)*2 )    (ones(numElemF,1)*3)      (zeros(numElemF,1))    (zeros(numElemF,1)) ...
                     %ElemNodes...
                     (1:(numElemF))'                         (2:numElemF+1)' ];
                     %MEBI truss elements
                   [ (ones(numElemT,1)*1 )    (ones(numElemT,1)*2)      (zeros(numElemT,1))     (zeros(numElemT,1)) ...
                     %ElemNodes...
-                    (numElemF + 1 : numElemF + numElemT)'	 	(numElemF + 2 : numElemF + numElemT + 1)'  ] 
-                ] ;  
+                    (numElemF + 1 : numElemF + numElemT)'	 	(numElemF + 2 : numElemF + numElemT + 1)'  ]
+                ] ;
 %md
 %mdFill conec cell with auxConecElem:
 for i =  1:numElemF + numElemT
-    mesh.conecCell{3 + i,1} = auxConecElem(i,:);
-end      
+    mesh.conecCell{2 + i,1} = auxConecElem(i,:);
+end
 %md
 %md### analysisSettings
 %md
-analysisSettings.deltaT        =   0.1 /2 ;
-analysisSettings.finalTime     =   4.13 ;
-% analysisSettings.methodName    = 'newmark' ;
-% analysisSettings.alphaNM       =   0.25 ;
-% analysisSettings.deltaNM       =   0.5  ;
+analysisSettings.deltaT        =   0.025 ;
+analysisSettings.finalTime     =   1.0 ;
 analysisSettings.methodName    = 'alphaHHT' ;
 analysisSettings.alphaHHT      =   -0.05   ;
 analysisSettings.stopTolDeltau =   1e-6 ;
@@ -129,4 +128,4 @@
 labx = xlabel('time (s)');   laby = ylabel('uT_z (m)') ;
 set(gca, 'linewidth', lw2, 'fontsize', plotfontsize )
 set(labx, 'FontSize', plotfontsize); set(laby, 'FontSize', plotfontsize) ;
-% print('output/rightAngleCantilever_uA_z.png','-dpng')
+% print('output/rightAngleCantilever_uA_z.png','-dpng')
diff --git a/makeAndPreviewDocs.sh b/makeAndPreviewDocs.sh
diff --git a/src/ONSAS_init.m b/src/ONSAS_init.m
@@ -19,7 +19,7 @@
 function [ modelCurrSol, modelProperties, BCsData ] = ONSAS_init( materials, elements, boundaryConds, initialConds, mesh, analysisSettings, otherParams )
 
 %md sets the current version
-ONSASversion = '0.2.3'  ;
+ONSASversion = '0.2.4'  ;
 
 %md set defaults
 [ materials, elements, boundaryConds, analysisSettings, otherParams ] = setDefaults( materials, elements, boundaryConds, analysisSettings, otherParams ) ;

diff --git a/src/vtk/frameVtkData.m b/src/vtk/frameVtkData.m
@@ -6,7 +6,7 @@
   vtkNodalDisps   = [] ;
   vtkNormalForces = [] ;
 
-  nPlotSubElements = 10 ; % number of plot subsegments
+  nPlotSubElements = 1 ; % number of plot subsegments
   counterNodes     = 0 ;
   nelem            = size(Conec,1) ;
 
@@ -82,21 +82,13 @@
       vtkConec             = [ vtkConec ;     Conecvtk ] ;
       vtkNodalDisps        = [ vtkNodalDisps; Dispsvtk ] ;
 
-      counterNodes = counterNodes + 8 ;
+      counterNodes = counterNodes + (size(Conecvtk,2)-1) ;
 
     end % for plot points
 
   end % for elements
 
 
-
-
-
-
-
-
-
-
   %
   %
   %       if j==1 % compute the first face of the first subelement

diff --git a/src/vtk/vtkBeam2SolidConverter.m b/src/vtk/vtkBeam2SolidConverter.m
@@ -21,38 +21,29 @@
   coordsElemNodes, dispsElem, coordLocSubElem, dispLocIni, dispLocEnd, ...
   locRotIni, locRotEnd, sectPar, Rr, R0 ) ;
 
-  typeSolid = sectPar(1) ;
+typeSolid = sectPar(1) ;
 
-  dispIniSection = dispsElem(1:2:5) ;    dispEndSection = dispsElem(7:2:11) ;
+dispIniSection = dispsElem(1:2:5) ;    dispEndSection = dispsElem(7:2:11) ;
 
-  if typeSolid == 12 % vtkHexa
+dispLocMed = 0.5 * ( dispLocIni + dispLocEnd ) ;
 
-     by = sectPar(2) ;  bz = sectPar(3) ;
+ex = [1 0 0] ;
+ey = [0 1 0] ;
+ez = [0 0 1] ;
 
-     Nodesvtk = zeros( 8,3 ) ;  Conecvtk = zeros( 8,1 ) ;
+if typeSolid == 12 % vtkHexa
 
-     % coordinates of nodes 1 and 2 in deformed configuration
-%5     node1Def = nodesCoords(1:3) + nodalDisp(1:3) ;
-%     node2Def = nodesCoords(4:6) + nodalDisp(4:6) ;
+   by = sectPar(2) ;  bz = sectPar(3) ;
 
-%     locglos = beamRefConfRotMat( (node2Def - node1Def ) ) ;
+   Nodesvtk = zeros( 8,3 ) ;  Conecvtk = zeros( 8,1 ) ;
 
-%     ex = locglos(:,1)' ;
-%     ey = locglos(:,2)' ;%
-%     ez = locglos(:,3)' ;
-     ex = [1 0 0] ;
-     ey = [0 1 0] ;
-     ez = [0 0 1] ;
-
-    % 1- compute the vectors of the section in Rr coords
-    % matrix with coords of four vertices of cross section to be plotted
-    matrixSectionIni = [ -ey*by*.5-ez*bz*.5 ; ...
-                         +ey*by*.5-ez*bz*.5 ; ...
-                         +ey*by*.5+ez*bz*.5 ; ...
-                         -ey*by*.5+ez*bz*.5 ] ;
-    matrixSectionEnd = matrixSectionIni ;
-    matrixSectionMed = [] ;
-  end
+  % 1- compute the vectors of the section in Rr coords
+  % matrix with coords of four vertices of cross section to be plotted
+  matrixSectionIni = [ -ey*by*.5-ez*bz*.5 ; ...
+                       +ey*by*.5-ez*bz*.5 ; ...
+                       +ey*by*.5+ez*bz*.5 ; ...
+                       -ey*by*.5+ez*bz*.5 ] ;
+  matrixSectionEnd = matrixSectionIni ;
 
   % 2- apply the local rotation for the initial and final sections
   matrixRotatedSectionIni = ( expon( locRotIni ) * matrixSectionIni' )' ;
@@ -78,65 +69,70 @@
   matrixRefIni = ( R0 * (matrixSectionIni + [ coordLocSubElem(1) 0 0] )')' + coordsElemNodes(1:3)' ;
   matrixRefEnd = ( R0 * (matrixSectionEnd + [ coordLocSubElem(2) 0 0] )')' + coordsElemNodes(1:3)' ;
 
-  NodesRefvtk = [ matrixRefIni ; matrixRefEnd ] ;
-
-  Dispsvtk = Nodesvtk - NodesRefvtk ;
-
+  NodesRefvtk = [ matrixRefIni; matrixRefEnd ] ;
 
-return
-	R = sectPar(2) / 2 ;
-  NodesDef  = nodesCoords + [ Ue(1:6:end) Ue(3:6:end) Ue(5:6:end) ] ;
-  rotsMat   = 			        [ Ue(2:6:end) Ue(4:6:end) Ue(6:6:end) ] ;
+elseif typeSolid == 25 % vtkQuadHexa
+  R = sectPar(2) / 2 ;
 
-  [~, locglos] = beamParameters( NodesDef ) ;
+  Nodesvtk = zeros( 20,3 ) ;  Conecvtk = zeros( 20,1 ) ;
 
+  matrixSectionIni = [ -ey*R/sqrt(2)-ez*R/sqrt(2) ; ...
+											 +ey*R/sqrt(2)-ez*R/sqrt(2) ; ...
+											 +ey*R/sqrt(2)+ez*R/sqrt(2) ; ...
+											 -ey*R/sqrt(2)+ez*R/sqrt(2) ] ;
+  matrixSectionEnd = matrixSectionIni ;
 
-  ex = locglos(:,1)' ;
-	ey = locglos(:,2)' ;
-	ez = locglos(:,3)' ;
+  matrixSectionCurvIni = [ 0-ez*R      ; ...
+                      +ey*R-0      ; ...
+                      0+ez*R       ; ...
+                      -ey*R+0      ] ;
 
-	%
-	matSecNodesExtVertices   = [ -ey*R/sqrt(2)-ez*R/sqrt(2) ; ...
-															 +ey*R/sqrt(2)-ez*R/sqrt(2) ; ...
-															 +ey*R/sqrt(2)+ez*R/sqrt(2) ; ...
-															 -ey*R/sqrt(2)+ez*R/sqrt(2) ] ;
+  matrixSectionCurvEnd = matrixSectionCurvIni ;
 
-	matSecNodesInter         = matSecNodesExtVertices ;
+  matrixSectionMed = matrixSectionIni ;
 
-	matSecNodesExtInter      = [             0-ez*R       ; ...
-                                       +ey*R-0           ; ...
-																				   0+ez*R       ; ...
-                                       -ey*R+0           ] ;
+  % 2- apply the local rotation for the initial and final sections
+  matrixRotatedSectionIni     = ( expon( locRotIni )                  * matrixSectionIni'     )' ;
+  matrixRotatedSectionEnd     = ( expon( locRotEnd )                  * matrixSectionEnd'     )' ;
+  matrixRotatedSectionCurvIni = ( expon( locRotIni )                  * matrixSectionCurvIni' )' ;
+  matrixRotatedSectionCurvEnd = ( expon( locRotEnd )                  * matrixSectionCurvEnd' )' ;
+  matrixRotatedSectionMed     = ( expon( (locRotIni+locRotEnd)*0.5  ) * matrixSectionMed'     )' ;
 
+  % 3- add local displacement and position
+  matrixDisplacedSectionIni     = matrixRotatedSectionIni     + [ coordLocSubElem(1) 0 0] + dispLocIni' ;
+  matrixDisplacedSectionEnd     = matrixRotatedSectionEnd     + [ coordLocSubElem(2) 0 0] + dispLocEnd' ;
+  matrixDisplacedSectionCurvIni = matrixRotatedSectionCurvIni + [ coordLocSubElem(1) 0 0] + dispLocIni' ;
+  matrixDisplacedSectionCurvEnd = matrixRotatedSectionCurvEnd + [ coordLocSubElem(2) 0 0] + dispLocEnd' ;
+  matrixDisplacedSectionMed     = matrixRotatedSectionMed     + [ (coordLocSubElem(1)+coordLocSubElem(2))*.5 0 0] + dispLocMed' ;
 
-	% Rot section 1
-	matSecNodesExtVerticesR  = ( expon( vecrotNode1 ) * matSecNodesExtVertices' )'  ;
-	matSecNodesExtInterR     = ( expon( vecrotNode1 ) * matSecNodesExtInter' )'     ;
+  % 4- apply Rr' change basis matrix
+  matrixRotatedSectionIni     = ( Rr * matrixDisplacedSectionIni'     )' ;
+  matrixRotatedSectionEnd     = ( Rr * matrixDisplacedSectionEnd'     )' ;
+  matrixRotatedSectionCurvIni = ( Rr * matrixDisplacedSectionCurvIni' )' ;
+  matrixRotatedSectionCurvEnd = ( Rr * matrixDisplacedSectionCurvEnd' )' ;
+  matrixRotatedSectionMed     = ( Rr * matrixDisplacedSectionMed'     )' ;
 
-	nodeExtVertices1  = ones(4,1) * NodesDef(1,:) + matSecNodesExtVerticesR    ;
-	nodeExtInter1     = ones(4,1) * NodesDef(1,:) + matSecNodesExtInterR       ;
+  defPosIniSec = coordsElemNodes(1:3) + dispIniSection ;
 
-	% Rot section 2
-	matSecNodesExtVerticesR  = ( expon( vecrotNode2 ) * matSecNodesExtVertices' )'  ;
-	matSecNodesExtInterR     = ( expon( vecrotNode2 ) * matSecNodesExtInter' )'     ;
+  % 5- add nodal displacements in e1,e2,e3 system
+  nodesVtkSectionIni     = matrixRotatedSectionIni     + defPosIniSec' ;
+  nodesVtkSectionEnd     = matrixRotatedSectionEnd     + defPosIniSec' ;
+  nodesVtkSectionCurvIni = matrixRotatedSectionCurvIni + defPosIniSec' ;
+  nodesVtkSectionCurvEnd = matrixRotatedSectionCurvEnd + defPosIniSec' ;
+  nodesVtkSectionMed     = matrixRotatedSectionMed     + defPosIniSec' ;
 
-	nodeExtVertices2  = ones(4,1) * NodesDef(2,:) + matSecNodesExtVerticesR     ;
-	nodeExtInter2     = ones(4,1) * NodesDef(2,:) + matSecNodesExtInterR        ;
+  Nodesvtk = [ nodesVtkSectionIni ;     nodesVtkSectionEnd; ...
+               nodesVtkSectionCurvIni;  nodesVtkSectionCurvEnd;  nodesVtkSectionMed ] ;
+  Conecvtk = [ 25 0:19 ] ; % in vtk indexation (from 0)
 
-  % Rot intermediate section
-	vecrot = ( vecrotNode1 + vecrotNode2 ) / 2            ;
-	matSecNodesInterR = ( expon( vecrot) * matSecNodesInter' )'               ;
-	nodeInt           = ones(4,1) * ( NodesDef(1,:) + NodesDef(2,:) ) / 2 + matSecNodesInterR ; % Interpolated linearly ...
+  matrixRefIni = ( R0 * (matrixSectionIni + [ coordLocSubElem(1) 0 0] )')' + coordsElemNodes(1:3)' ;
+  matrixRefEnd = ( R0 * (matrixSectionEnd + [ coordLocSubElem(2) 0 0] )')' + coordsElemNodes(1:3)' ;
+  matrixRefCurvIni = ( R0 * (matrixSectionCurvIni + [ coordLocSubElem(1) 0 0] )')' + coordsElemNodes(1:3)' ;
+  matrixRefCurvEnd = ( R0 * (matrixSectionCurvEnd + [ coordLocSubElem(2) 0 0] )')' + coordsElemNodes(1:3)' ;
+  matrixRefMed = ( R0 * (matrixSectionMed + [ (coordLocSubElem(1)+coordLocSubElem(2))*.5 0 0] )')' + coordsElemNodes(1:3)' ;
 
-  % Nodes
-	Nodesvtk = [ Nodesvtk         ; ...
-               nodeExtVertices1 ; ...
-               nodeExtVertices2 ; ...
-               nodeExtInter1    ; ...
-               nodeExtInter2    ; ...
-               nodeInt          ] ;
+  NodesRefvtk = [ matrixRefIni ;  matrixRefEnd; matrixRefCurvIni; matrixRefCurvEnd; matrixRefMed ] ;
 
-	% Connectivity
-	Conecvtk = [ Conecvtk ; 25 1:20 ] ;
+end
 
-%end
+Dispsvtk = Nodesvtk - NodesRefvtk ;
diff --git a/src/vtk/vtkDataConversion.m b/src/vtk/vtkDataConversion.m
@@ -52,6 +52,7 @@
         = trussVtkData( modP.Nodes, modP.Conec( elemIndsElemType, 5:end ), ...
         elemTypeGeom, modS.U ) ;
 
+
     elseif strcmp( elemTypeString, 'frame' )
 
       [ currVtkNodes, currVtkConec, currVtkNodalDisps, vtkNormalForces ] ...
@@ -87,18 +88,25 @@
 
             % if nargout > numminout
               % add the tetrahedron vtk cell type to the first column
-      currVtkConec = [ 10*ones( nelems, 1 )     modP.Conec(:, 5:8 )-1 ] ;
-            elem2VTKCellMap = (1:nelems)' ; % default: i-to-i . Columns should be changed.
-            % end
-
+      currVtkConec    = [ 10*ones( nelems, 1 )     modP.Conec(:, 5:8 )-1 ] ;
+      elem2VTKCellMap = (1:nelems)' ; % default: i-to-i . Columns should be changed.
     end % if: type
 
     % add entries from current element type
     vtkNodes      = [ vtkNodes ;  currVtkNodes ]           ;
-    vtkConec      = [ vtkConec ;  [currVtkConec(:,1) currVtkConec(:,2:end)+totalNodes]]           ;
+
+    if size(vtkConec,1) > 0
+      vtkConec( ...
+        (size(vtkConec,1)+1):(size(vtkConec,1)+size(currVtkConec,1)), 1:(size(currVtkConec,2)) ) ...
+        = [currVtkConec(:,1) currVtkConec(:,2:end)+totalNodes] ;
+    else
+      vtkConec = [currVtkConec(:,1) currVtkConec(:,2:end)+totalNodes] ;
+    end
+
     vtkNodalDisps = [ vtkNodalDisps ;  currVtkNodalDisps ] ;
 
-		totalNodes = totalNodes + size(currVtkNodes,1) ; 
+		totalNodes = totalNodes + size(currVtkNodes, 1) ;
+
 
   end % for: elemTypeInds