APL(NARS), 58 chars

{(g m a)←⍎¨⍵⊂⍨⍵≠'.'⋄8∣⌊0.5+3.6÷⍨30∣(11×19∣a-m≤2)+g+12∣m-3}

The function it is based of one Internet site speak how the ancients calculate moon phases and I change something for golfing. The site i have seen says +- what the site of the link says

https://farelaboratorio.accademiadellescienze.it/esperimenti/scienze/65

with Epatta as that.

Ancients consider moon calendar begin to 3th month (m≤2 m-3), "Marzo" March, while 1th month Gen and 2th month Feb are with the year that is before. Than someone observed that every 19 years the phases repeated with the same days. Number 11 i don't know why there is in the formula.

It is possible that i make wrong something and the error result more for some date. The site claim formula would be ok for phases, from 1799..2199 but I have had changed it and would result ok from 1900 to 2199.

From here Translated by Google and rewrote for changes.

I took the lunar month composed of about 29.53 days round to 30 days, I divided it by 8 equal parts, so I got 8 sub-intervals numbered from 0 to 7, they are the 8 lunar phases (if I understood correctly).

So the day represented by day/month/year will also be also a day of the lunar month represented by a number N between 0 and 29, which in turn can be represented by another number that expresses the phase of the moon that goes from 0 to 7 depending on whether the number N falls in one of the 8 equal parts in which the lunar month 0..29 is divided from part 0 to part 7 (round toward the number of phase more near)

The function above has as input the "day.month.year" as one string, and as output the value of the phase from 0 to 7 above defined. Phases and their names:

   0       1             2                3             4           5                6             7 
('New')('Waxing')('First Quarter')('Waxing Gibbous')('Full')('Waning Gibbous')('Last Quarter')('Waning') 

Test month of January 2016 days from 7 to 16, month of February 2016 from 1 to 11, and other days:

   H←{(g m a)←⍎¨⍵⊂⍨⍵≠'.'⋄8∣⌊0.5+3.6÷⍨30∣(11×19∣a-m≤2)+g+12∣m-3}
  ⍪{⍵,H ⍵}¨("14.12.2016")("16.10.1983")("4.7.1976")("28.11.1970")("7.1.1970")("27.4.2025")("1.1.1900")("31.1.1900")
14.12.2016 4 
16.10.1983 3 
4.7.1976 2   
28.11.1970 0 
7.1.1970 0   
27.4.2025 0  
1.1.1900 0   
31.1.1900 0  
  {⍵, H (⍕⍵),".1.2016"}¨7..16 
7 0  8 0  9 0  10 0  11 1  12 1  13 1  14 1  15 2  16 2 
  {⍵, H (⍕⍵),".2.2016"}¨1..11 
1 6  2 7  3 7  4 7  5 0  6 0  7 0  8 0  9 0  10 1  11 1 

This below function K would be ok from 1500 to 3099 and use one table found in the site link above, it has input string "day.month.year", it would return one number that represent the percentage of the illumination of the moon, and one other number represent the phase.

  K←{(g m a)←⍎¨⍵⊂⍨⍵≠'.'⋄a-←m≤2⋄v←1,1,0,0,-1,1,1,2,3,2,3,4,4,4,5,5⋄(⌊{k≤15:k⋄30-k}×20÷3),8∣⌊0.5+3.625÷⍨k←30∣(v[¯14+⌊a÷100]+11×19∣a)+g+1+12∣m-3}
  ⍪{⍵,K ⍵}¨("14.12.2016")("16.10.1983")("4.7.1976")("28.11.1970")("7.1.1970")("27.4.2025")("1.1.1900")("31.1.1900")
14.12.2016 100 4 
16.10.1983 66 3  
4.7.1976 53 2    
28.11.1970 6 0   
7.1.1970 6 0     
27.4.2025 6 0    
1.1.1900 0 0     
31.1.1900 0 0    

A number equal to or close to 100% represents a full moon, while a number equal to or close to 0% represents a new moon, a number equal to or close to 50% represents a quarter moon, first quarter moon if it is around phase 2, last quarter moon if it is around phase 6.

I hope that is as below the speach for errors:

The error for H for date in years 1900..2199 should be 0 for phase (as here defined), +-4 days for the day has in their time the hour when there is 100% moon full, or 50% moon full, or 0% moon full, +-1 the phase if it is calculated as it would be 1 only day, example "phase full moon" as only the day there is full moon.

The error for K for date in years 1500..3099 should be 0 for phase (as here defined), +-1 day for the day has in their time the hour when there is 100% moon full, or 50% moon full, or 0% moon full, because the day will be the day with the right percentage of the day, or the one near to it in the left, or the one near to it in the right. Example for 10-1-2016 and 8-2-2016 New Moon:

  {⍵, K (⍕⍵),".1.2016"}¨7..16 
7 13 0  8 6 0  9 0 0  10 6 0  11 13 1  12 20 1  13 26 1  14 33 1  15 40 2  16 46 2 
  {⍵, K (⍕⍵),".2.2016"}¨1..11 
1 46 6  2 40 7  3 33 7  4 26 7  5 20 7  6 13 0  7 6 0  8 0 0  9 6 0  10 13 1  11 20 1

10-1-2016 is near in the right to 9-1-2016 (the day of precise percentage) so the error 1 day, and 8-2-2016 is the value with right percentage so no error.

Python 2 3, 255 204 180 178 bytes

This is answer is inaccurate by a day or two in several places, including for some of the test cases, though I was told that some inaccuracy was acceptable. At any rate, the motion of the moon is never very exact and this function remains generally correct (or at least, it doesn't vary too far).

Edit: In the course of correcting my code and making it more accurate, I have golfed it down considerably.

Edit: This code is now a one-line Python 3 program. (Credit to TimmyD for the name "magic numbers")

p,q,r=(int(i)for i in input().split("-"));t=q<3;y=p-2000-t;i,j=divmod(((r+(153*(q+12*t-3)+2)//5+365*y+y//4-y//100+y//400+11010)*86400-74100)%2551443,637861);print((86400<=j)+2*i)

Ungolfed:

def jul(p,q,r):
    '''
    The Julian Day Number (JDN) of the input minus the JDN of January 7, 1970,
    the first new moon after January 1, 1970.
    '''
    t=q<3
    y=p-2000-t  # -2000 years to push day 0 to January 1, 2000
    return r+(153*(q+12*t-3)+2)//5+365*y+y//4-y//100+y//400+11010
    # +11010 days to push day 0 to January 7, 1970

def moon(s):
    '''
    Input format: yyyy-mm-dd

    An attempt at explaining the "magic numbers"
    - 29.53059 days is close to 2551443 seconds, so I used that
    - The offset of +12300 seconds because the new moon of 1970-01-07 was at 2035 UTC 
      or 12300 seconds before midnight. For those of you saying that this pushes 
      the beginning of my calendar to 2035, *6* January 1970, yes it does.
      But if I need to the calendar to recognize 1970-01-07 as the day of the new moon 
      which means that midnight needed to be a positive number of seconds, 0 <= x < 86400.
      Basically, I hacked it together, and +12300 worked.        
    '''
    d = 86400
    p,q,r = map(int, s.split("-"))
    z=(jul(p,q,r)*d+12300)%2551443  # 2551443 is about the number of seconds in a lunar month
    div, mod = divmod(z, 637861)    # 637861 seconds is about a quarter of a lunar month
                                    # div is what part of the lunar month this is (0 - 3)
                                    # mod is seconds since the start of the main phase
    return 2*div + (86400 <= mod)   # 2*div for the main phase, and 
                                    # is mod >= the number seconds in a day?
                                    # (+0 if within a day of the main phase, +1 otherwise)