Merge pull request #3 from jmpearl/Update-Integrator

Update integrator

scripts/Reproducibility/masconPaper/Fig04_LayeredMascons.m

@@ -0,0 +1,63 @@
+clear all; close all; clc;
+addpath(genpath('/home/jmpearl/GravityModels/CLEO'))
+
+Mu = 1.0;
+
+meshfile = 'Eros_46906.obj'; % mesh file to load
+sm = SurfaceMesh(meshfile);
+sm.coarsen(15000);
+
+vm = VolumeMesh(sm);
+vm.initializeFromSurfaceIteration(1/3)
+
+masconModel = MasconModel();
+masconModel.initializeFromVolumeMesh(vm,Mu,"excludesurface")
+
+vm.addCellField(vm.cellVolumes(),'cellVolumes')
+vm.writeVTK("Fig4_volumeMesh.vtk")
+
+
+% Fig 10 (entry for one body)
+%-----------------------------
+FS=16; LW=1; MS=8;
+thresh = 2000;
+markerScaleFactor = 24;
+figure(4)
+hold on
+
+% plot faces in certain range x
+c = sm.faceCentroids();
+for i = 1:sm.numFaces
+ if (c(i,3) < thresh && c(i,3) > -thresh)
+ pts = [sm.coordinates(sm.faces(i,:),:);...
+ sm.coordinates(sm.faces(i,1),:)];
+ plot3(pts(:,1),pts(:,2),pts(:,3),'k-')
+ end
+end
+
+c = masconModel.coordinates;
+mu = masconModel.mu;
+
+resolution = mu.^(1/3);
+MS = resolution./max(resolution)*markerScaleFactor;
+color = 1.0-(MS-min(MS))/max(max(MS)-min(MS),mean(MS)/10000);
+colorVecMax = [1,1,0];
+colorVecMin = [0,0,1];
+for i = 1:masconModel.numElements
+ if (c(i,3) < thresh && c(i,3) > -thresh)
+ colorVeci = (1-color(i))*colorVecMax + color(i)*colorVecMin;
+ MSi = MS(i);
+ plot3(masconModel.coordinates(i,1),...
+ masconModel.coordinates(i,2),...
+ masconModel.coordinates(i,3),'ko','MarkerFaceColor',colorVeci,'MarkerSize',MSi)
+ end
+end
+set(gcf,'Color',[1,1,1]);
+set(gca,'FontSize',FS)
+set(gca,'TickLabelInterpreter','latex')
+daspect([1,1,1])
+view([0,0,1])
+box off
+axis off
+
+%[f1, cdata] = myaa([4, 2]); imwrite(cdata, [fileName], 'png');

src/GravityModels/ApproximatePolyhedralModel.m

@@ -79,7 +79,7 @@ function initializeFromMesh(obj,mesh,Mu,quadratureRule)
  for i = 1:size(p,1)
 
  rhat = obj.positions-p(i,:); 
- rhat = rhat./sqrt(rhat(:,1).^2+rhat(:,2).^2+rhat(:,3).^2);
+ rhat = rhat./max(sqrt(rhat(:,1).^2+rhat(:,2).^2+rhat(:,3).^2),0.2*obj.res(i));
  %rhat(isnan(rhat)|isinf(rhat))=0;
  potential(i,1) = -1/2 * rhat(:)' * obj.GrhoAn(:); 
 

src/GravityModels/MasconModel.m

@@ -194,6 +194,19 @@ function initializeExtendedTetrahedra(obj, mesh, Mu, numLayers, method)
  % get our resolution
  obj.calculateResolution(volumes);
  end
+
+ function c = centroid(obj)
+ % moves distributions so centroid is at origin
+ %------------------------------------------------------------------
+ firstMoment = sum(obj.mu,1);
+ secondMoment = sum(obj.mu.*obj.coordinates,1);
+ c = secondMoment/firstMoment;
+ end
+ function center(obj)
+ % moves distributions so centroid is at origin
+ %------------------------------------------------------------------
+ obj.coordinates=obj.coordinates-obj.centroid();
+ end
 
  function calculateResolution(obj,volumes)
  obj.res = 0.5*volumes.^(1/3);

src/Integrator/Integrator.m

@@ -6,11 +6,14 @@
  gravityModel; % models applied in BFF frame
  thirdBodyModel; % interial frame 3rd body acceleration
  SRPModel; % solar radiation pressure model
- omega;  % rotation rate (rad/sec)
+ omega=0.0; % rotation rate (rad/sec)
 
  integrator; % function handle for integrator
  mode = 0; % 0-inertial, frame 1-BFF, 2-variational Eqns
  odeOptions; % options from odeset
+ inertialThirdBodyFrame=true % option to swtich third body 
+ % acceleration to be BFF integrated to
+ % represent tidal locked bodies
  end
 
  methods
@@ -101,7 +104,11 @@ function setOdeOptions(obj,options)
  a = obj.gravityModel.acceleration(P_bff)*Rot'; 
 
  if ~isempty(obj.thirdBodyModel)
- a = a + obj.thirdBodyModel.acceleration(P_inertial);
+ if obj.inertialThirdBodyFrame
+ a = a + obj.thirdBodyModel.acceleration(P_inertial);
+ else
+ a = a + obj.thirdBodyModel.acceleration(P_bff)*Rot';
+ end
  end
  if ~isempty(obj.SRPModel)
  a = a + obj.SRPModel.acceleration(P_inertial);

src/Mesh/SurfaceMesh.m

@@ -1565,10 +1565,10 @@ function setResolution(obj,newResolution)
 
  newNumFaces = obj.surfaceArea/newResolution^2;
 
- if newNuFaces < 20
+ if newNumFaces < 20
  warning('very low face count might fail')
  end
-  warni
+
  obj.setNumFaces(newNumFaces);
 
  end
@@ -3028,7 +3028,7 @@ function checkValid(obj)
  elseif size(obj.faceFields{i}.data,2)==3
  vtkDataHeader = ['VECTORS ',obj.faceFields{i}.name,' float\n'];
  fprintf(fileID,vtkDataHeader);
- fprintf(fileID,formatSpecVectorFloat,obj.faceFields{i}.data);
+ fprintf(fileID,formatSpecVectorFloat,obj.faceFields{i}.data');
  else
  warning(['skipping faceData{',num2str(i),'} not scalar or vector'])
  end
@@ -3054,7 +3054,7 @@ function checkValid(obj)
  elseif size(obj.nodeFields{i}.data,2)==3
  vtkDataHeader = ['VECTORS ',obj.nodeFields{i}.name,' float\n'];
  fprintf(fileID,vtkDataHeader);
- fprintf(fileID,formatSpecVectorFloat,obj.nodeFields{i}.data);
+ fprintf(fileID,formatSpecVectorFloat,obj.nodeFields{i}.data');
  else
  warning(['skipping faceData{',num2str(i),'} not scalar or vector'])
  end

src/Mesh/VolumeMesh.m

@@ -292,6 +292,7 @@ function initializeFromSurfaceIteration(obj,ratio,numVerticesInitial)
 
  % offset the surface mesh inward coarsening along the way
  while tempMesh.volume > 0.1*obj.surfaceMesh.volume
+ %tempMesh.addFaceField(tempMesh.faceAreas(),'faceAreas')
  %tempMesh.writeVTK(['tempMesh',num2str(i),'.vtk'])
  i=i+1;
  % reduce vertex count 
@@ -845,7 +846,51 @@ function smooth(obj,numIterations)
  boolOutput(obj.numBoundaryVertices+1:end,:)=0;
 
  end
-
+
+ function [edgeLengths] = cellEdgeLengths(obj)
+ % edge lengths for each cell
+ %------------------------------------------------------------------
+ % Outputs:
+ % centroids - coordinates of centroids
+ % volumes --- volumes of each cell
+ %------------------------------------------------------------------
+ p1 = obj.coordinates(obj.cells(:,1),:);
+ p2 = obj.coordinates(obj.cells(:,2),:); 
+ p3 = obj.coordinates(obj.cells(:,3),:); 
+ p4 = obj.coordinates(obj.cells(:,4),:); 
+ v1 = vecnorm(p2-p1,2,2);
+ v2 = vecnorm(p3-p1,2,2);
+ v3 = vecnorm(p4-p1,2,2);
+ v4 = vecnorm(p2-p4,2,2);
+ v5 = vecnorm(p3-p2,2,2);
+ v6 = vecnorm(p4-p3,2,2);
+
+ edgeLengths = [v1,v2,v3,v4,v5,v6];
+ end
+ function [faceAreas] = cellFaceAreas(obj)
+ % areas of the foue faces
+ %------------------------------------------------------------------
+ % Outputs:
+ % centroids - coordinates of centroids
+ % volumes --- volumes of each cell
+ %------------------------------------------------------------------
+ p1 = obj.coordinates(obj.cells(:,1),:);
+ p2 = obj.coordinates(obj.cells(:,2),:); 
+ p3 = obj.coordinates(obj.cells(:,3),:); 
+ p4 = obj.coordinates(obj.cells(:,4),:); 
+ v1 = p2-p1;
+ v2 = p3-p1;
+ v3 = p4-p1;
+ v4 = p2-p4;
+ v5 = p3-p2;
+ %v6 = p4-p3;
+
+ A1 = 0.5*vecnorm(cross(v1,v2),2,2);
+ A2 = 0.5*vecnorm(cross(v1,v3),2,2);
+ A3 = 0.5*vecnorm(cross(v2,v3),2,2);
+ A4 = 0.5*vecnorm(cross(v4,v5),2,2);
+ faceAreas = [A1,A2,A3,A4];
+ end
  function [centroids,volumes] = cellCentroids(obj)
  % centroids of cells
  %------------------------------------------------------------------
@@ -1614,6 +1659,23 @@ function writeVTK(obj,fileName,precision)
 
  end
 
+ function plotVertexCells(obj,i)
+
+ figure()
+ hold on
+ for j = 1:obj.numCells
+ if any(obj.cells(j,:)==i)
+ p1 = obj.coordinates(obj.cells(j,1),:);
+ p2 = obj.coordinates(obj.cells(j,2),:); 
+ p3 = obj.coordinates(obj.cells(j,3),:); 
+ p4 = obj.coordinates(obj.cells(j,4),:); 
+ X = [p1;p2;p3;p4;p1;p3;p4;p2];
+ plot3(X(:,1),X(:,2),X(:,3),'k-')
+ %break
+ end
+ end
+ daspect([1,1,1])
+ end
  function plotCellNodes(obj,i)
 
  figure(1)