import importDAE, Part, ObjectsFem, Fem, FemGmshTools, FemGui, FemToolsCcx

PRESSURE_MODE = 1
MAX_DEPTH = 8.0
HEIGHT = 6.48
GROUND_DENSITY = 1500.0
GMSH_RES_METERS = 2.0
PRESSURE_LAYERS = 100

importDAE.open(u"/tmp/tmp.dae")
doc = FreeCAD.getDocument("tmp")
gui_doc = Gui.getDocument("tmp")
mesh = doc.getObject("Mesh")
ground = doc.getObject("Mesh001")
Gui.SendMsgToActiveView("ViewFit")

print "Object loaded.  Making shape..."
part = doc.addObject("Part::Feature", "Hús part")
shape_tmp = Part.Shape()
shape_tmp.makeShapeFromMesh(mesh.Mesh.Topology, 0.100000)
part.Shape = shape_tmp
part.purgeTouched()
delete shape_tmp

# solidify
print "Solidifying..."
solid_tmp = part.Shape
solid_tmp = Part.Solid(solid_tmp)
solid_part = doc.addObject("Part::Feature", "Hús solid part")
solid_part.Label = "Hús solid part"
solid_part.Shape = solid_tmp
delete solid_tmp

# refining
print "Refining..."
refined_part = doc.addObject("Part::Feature", "Hús refined part")
refined_part.Shape = solid_part.Shape.removeSplitter()

# analysis
analysis_object = ObjectsFem.makeAnalysis("Analysis")
FemGui.setActiveAnalysis(doc.Analysis)

# gmsh
print "Computing gmsh..."
gmsh = ObjectsFem.makeMeshGmsh('Hús FEM')
gmsh.Part = refined_part
gmsh.CharacteristicLengthMax = 1000.0 * GMSH_RES_METERS
gmsh.CharacteristicLengthMin = 0.0
gmsh_mesh = FemGmshTools.FemGmshTools(gmsh)
error = gmsh_mesh.create_mesh()
print error
doc.Analysis.Member += [gmsh]
doc.recompute()

# solver
print "Creating solver..."
solver_object = ObjectsFem.makeSolverCalculix("CalculiX")
solver_object.GeometricalNonlinearity = 'linear'
solver_object.ThermoMechSteadyState = True
solver_object.MatrixSolverType = 'default'
solver_object.IterationsControlParameterTimeUse = False
doc.Analysis.Member = doc.Analysis.Member + [solver_object]

# material
print "Adding material..."
material_object = ObjectsFem.makeMaterialSolid("SolidMaterial")
mat = material_object.Material
mat['Name'] = "Pumicecrete"
mat['YoungsModulus'] = "5000 MPa"
mat['PoissonRatio'] = "0.17"
mat['Density'] = "1000 kg/m^3"
material_object.Material = mat
doc.Analysis.Member = doc.Analysis.Member + [material_object]

# bounding
print "Bounding..."
MIN_X =  9999999999999.9
MAX_X = -9999999999999.9
MIN_Y =  9999999999999.9
MAX_Y = -9999999999999.9
MIN_Z =  9999999999999.9
MAX_Z = -9999999999999.9
for vertex in refined_part.Shape.Vertexes:
  if vertex.X > MAX_X:
    MAX_X = vertex.X
  if vertex.X < MIN_X:
    MIN_X = vertex.X
  if vertex.Y > MAX_Y:
    MAX_Y = vertex.Y
  if vertex.Y < MIN_Y:
    MIN_Y = vertex.Y
  if vertex.Z > MAX_Z:
    MAX_Z = vertex.Z
  if vertex.Z < MIN_Z:
    MIN_Z = vertex.Z

SCALAR_X = MAX_X - MIN_X
SCALAR_Y = MAX_Y - MIN_Y
SCALAR_Z = (MAX_Z - MIN_Z) / HEIGHT


# fixed_constraint
print "Adding anchors..."
fixed_constraint = ObjectsFem.makeConstraintFixed("Anchors")
keys = []
for face_id in range(len(refined_part.Shape.Faces)):
  face = refined_part.Shape.Faces[face_id]
  z = (face.CenterOfMass.z - MIN_Z) / SCALAR_Z
  if z < 2.7 and face.normalAt(0, 0)[2] < -0.9:
    keys.append( (refined_part, "Face%d" % (face_id + 1)) )

fixed_constraint.References = keys
doc.Analysis.Member = doc.Analysis.Member + [fixed_constraint]

def sign(p0, p1, p2):
  return (p0[0] - p2[0]) * (p1[1] - p2[1]) - (p1[0] - p2[0]) * (p0[1] - p2[1])

def point_in_triangle(p0, p1, p2, x, y):
  pt = (x, y)
  b1 = sign(pt, p0, p1) < 0.0
  b2 = sign(pt, p1, p2) < 0.0
  b3 = sign(pt, p2, p0) < 0.0
  return (b1 == b2) & (b2 == b3)

def distance(x0, y0, x1, y1):
  return ((x0-x1)**2 + (y0-y1)**2)**0.5

def ground_height(x, y):
  for facet in ground.Mesh.Facets:
    p = []
    for i in facet.Points:
      p.append( ((i[0] - MIN_X) / SCALAR_X, (i[1] - MIN_Y) / SCALAR_Y, (i[2] - MIN_Z) / SCALAR_Z) )
    if point_in_triangle(p[0], p[1], p[2], x, y):
      break
  d0 = distance(x, y, p[0][0], p[0][1])
  d1 = distance(x, y, p[1][0], p[1][1])
  d2 = distance(x, y, p[2][0], p[2][1])
  sum = (d0 + d1 + d2)
  z = (d0 * p[0][2] + d1 * p[1][2] + d2 * p[2][2]) / sum
#  print x, y, " : ", p[0][2], p[1][2], p[2][2], " : ", d0, d1, d2, " : ", z
  return z

print "Adding pressures..."

if PRESSURE_MODE:
  # pressure_constraint
  keys = []
  for i in range(PRESSURE_LAYERS):
    keys.append([])
  for face_id in range(len(refined_part.Shape.Faces)):
    face = refined_part.Shape.Faces[face_id]
    average_depth = 0.0
    for vertex in face.Vertexes:
      x = (vertex.X - MIN_X) / SCALAR_X
      y = (vertex.Y - MIN_Y) / SCALAR_Y
      z = (vertex.Z - MIN_Z) / SCALAR_Z
      ground_z = ground_height(x, y)
      depth = ground_z - z
#      print ground_z, z, depth
      if depth > 0:
        average_depth += depth
    average_depth /= len(face.Vertexes)
    if average_depth > 0.0:
#      print "Adding key: %d %f to %d" % (face_id + 1, average_depth, int(round(average_depth / MAX_DEPTH * len(keys))))
      keys[int(round(average_depth / MAX_DEPTH * len(keys)))].append( (refined_part, "Face%d" % (face_id + 1)) )

  for key_id in range(len(keys)):
    if len(keys[key_id]) == 0:
      continue
    depth = key_id * MAX_DEPTH / len(keys)
#    print "%d) Adding pressures for %f: " % (key_id, depth), keys[key_id]
    pressure = (9.81 * GROUND_DENSITY * depth) / 1000000.0
    pressure_constraint = ObjectsFem.makeConstraintPressure("Pressure %0.4f" % pressure)
    pressure_constraint.References = keys[key_id]
    pressure_constraint.Pressure = pressure
    pressure_constraint.Reversed = False
    doc.Analysis.Member = doc.Analysis.Member + [pressure_constraint]
    #break

else:
  # ATH: ekki tilbúið!
  for key_id in range(len(keys)):
    if len(keys[key_id]) == 0:
      continue
    depth = key_id * MAX_DEPTH / len(keys)
    force_per_m2 = (9.81 * GROUND_DENSITY * depth)
    m2_on_z = 0# ATH
    force = force_per_m2 / m2_on_z
    force_constraint = ObjectsFem.makeConstraintForce("Force %0.0f" % force)
    force_constraint.References = keys[key_id]
    force_constraint.Force = force
    force_constraint.Reversed = False
    force_constraint.Direction = (refined_part, [nafn]) # ATH
    doc.Analysis.Member = doc.Analysis.Member + [force_constraint]

print "Pressures added."

## View object
#femmesh_obj = doc.addObject('Fem::FemMeshObject', 'FEM view')
#femmesh_obj.FemMesh = gmsh
#doc.Analysis.Member = doc.Analysis.Member + [femmesh_obj]

# run the analysis
FemGui.setActiveAnalysis(doc.Analysis)
fea = FemToolsCcx.FemToolsCcx()
fea.update_objects()
message = fea.check_prerequisites()
if not message:
  fea.reset_all()
  fea.run()
  fea.load_results()
else:
  FreeCAD.Console.PrintError("Houston, we have a problem! {}\n".format(message))  # in report view
  print("Houston, we have a problem! {}\n".format(message))  # in python console

# run the analysis
for m in analysis_object.Member:
  if m.isDerivedFrom('Fem::FemResultObject'):
    result_object = m
    break

gui_doc.getObject("Mesh001").Visibility = False
gui_doc.getObject("Mesh").Visibility = False
gui_doc.getObject("Hús part").Visibility = False
gui_doc.getObject("Hús solid part").Visibility = False
gui_doc.getObject("Hús refined part").Visibility = False

gmsh.ViewObject.setNodeDisplacementByVectors(result_object.NodeNumbers, result_object.DisplacementVectors)
gmsh.ViewObject.applyDisplacement(10)